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Minimum time accelerations of aircraft turbofan engines are presented. 
The calculation of these accelerations is made by using a piecewise-linear 
engxne model, and a new algorithm based on nonlinear programming. Use of 
this model and algorithm allows such trajectories to be readily calculated 
on a digital computer with a minimal expenditure of computer time. 

The new algorithm may be used for solution of optimal control problems 
which are nonlinear in the state variables, and linear in the control vari- 
ables. Specifically, the most general problem considered is to 'miaimize a 
performance index subject to satisfaction of the system dynamic equations, a 
set of terminal constraints, and path inequality constraints. The perform- 
ance index, system equations, and path constraints are all linear in the con- 
trol variables. 

It^ is shown that the optimal control for such problems is bang-bang, 
except for possible singular arcs, which are not considered. The algorithm 
requires that a nominal bang-bang solution be found which satisfies the sys- 
tem dynamic equations and terminal constraints. Once such a feasible solu- 
tion has been found, influence functions are generated which determine if the 
necessary conditions for optimality have be‘£‘n satisfied. If not, additional 
control switches are needed. Nonlinear optimization (gradient search) tech- 
niques are then used to vary the control switching times in order to improve 
the solution. 

The algorithm is used to find minimum time acceleration histories for 
the FI 00 engine, a two-spool turbofan engine which powers the F15 and F16 


aircraft o A piecewise-linear engine model is used. The linearized model 
used at a given time in the trajectory is determined by calculating a normal- 
ised "distance” from the cxirrent state to the state at each of several equi- 
librium points; the model linearized about the "closest” equilibrium point 
is then used. Minimum time solutions are obtained, and the resulting con- 
trol histories are used as inputs to a nonlinear simulation of the FlOO engine 
to verify the accuracy of the piecewise linear solution. 

In addition to the transient results , the linear models are also used 
to find the control settings which minimize steady-state specific fuel con- 
sumption. 
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ABSTRACT 


Minimuin timfi accelerations of aircraft turbofan engines are pre- 
sented. The calculation of these accelerations is made by using a 
piecewise-linear engine model, and a new algorithm based on nonlinear 
programming. Use of this model and algorithm allows such trajectories 
to be readily calculated on a digital computer with a minimal expendi- 
ture of computer time. 

The new algorithm may be used for solution of optimal control prob- 
lems which are nonlinear in the state variables , and linear in; the con- 
trol variables. Specifically, the most general problem considered is to 
minimize a performance index subject to satisfaction of the system 
dynamic equations, a set of terminal constraints, and path inequality 
constraints. The performance index, system equations, and path con- 
straints are all linear in the control variables. 

It is shown that the optimal control for such problems is bang- 
bang, except for possible singular arcs , which are not considered. The 
algorithm requires that a nominal bang-bang solution be found which 
satisfies the system dynamic equations and terminal constraints. Once 
such a feasible solution has been found, influence functions are gener- 
ated which determine if the necessary conditions for optimality have 
been satisfied. If not , additional control switches are needed. Non- 
linear optiinization (gradient search) techniques are then used to vary 
the control switching times in order to improve the solution. 
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The algorithm is used to find minimum time acceleration his- 
tories for the FlOO engine, a two-spool turbofan engine which powers 
the F15 and F16 aircraft. A pie.cewise-linear engine model is used. 

The linearized model used at a given time in the trajectory is deter- 
mined by calculating a normalized "distance" from the current state 
to the state at each of several equilibrium points; the model linear- 
ized about the "closest" equilibrium point is then used. Minimum time 
solutions are obtained, and the resulting control histories are used 
as inputs to a nonlinear simulation of the FlOO engine to verify the 
accuracy of the piecewise linear solution. 

In addition to the transient results, the linear models are also 
used to find the control settings which minimize steady-state specifxc 
fuel consumption. 
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CHAPTER I . INTRODUCTION 
A. Motivation 

Modern, high-performance turbojet and turbofan engines are gen- 
erally equipped with one or more variable geometry features in order to 
provide maximum propulsive efficiency over a range of engine power set- 
tings and flight conditions. For example, the J-85 engine (a one-spool 
turbojet used in the F5 aircraft) has variable inlet guide vanes, and 
variable bleeds in three stages of its eight-stage compressor. The 
TF30 engine (a two-spool turbofan, used in the Fill and F14 aircraft) 
has variable bleeds in the low and high compressors. The FlOO engine 
(a two -spool turbo fan, used in the F15 and F16 aircraft) has variable 
fan inlet guide vanes and variable compressor stator vanes. Ejach pf 
these engines also has a variable-area exhaust nozzle and an after- 
burner. Variable area turbines, although not yet in operational use, 
have been tested on technology demonstrator engines. 

Propulsive efficiency is probably the most important measure of an 
aircraft engine’s performance. However, another important measure is 
the time required to accelerate from one thrust level to a higher 
thrust level. Engine acceleration is one of the functions of the 
engine control system, and may be accomplished via open-loop scheduling 
or closed-loop control. For each of the three engines referred to 
above, engine accelerations are accomplished by controlling fuel flow. 
The variable geometry features are kept on their steady-state schedules 
during the acceleration. 
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In this report, minimum-time acceleration histories are computed 
for the FlOO turbofan engine. Four control variables, i.e., fuel flow, 
exhaust nozzle area, inlet guide vane position, and compressor stator 
vane position, are utilized. 

B . Related Work 

In recent years, linear-quadratic regulator theory has been devel- 
oped for the design of multi-input, multi-output control systems. An 
account of the theory and application is given for example in refer- 
ence 1. Use of the theory has been facilitated by computer programs 
such as those described in references 2 and 3, which rapidly and effi- 
ciently calculate the optimal feedback control gains , given the system 
description and performance index. This theory has been applied re- 
cently to the design of control systems for aircraft gas turbine en- 
gines. In addition to the design of regulators, the problem of mini- 
mizing acceleration time has also been considered. 

Michael and Farrar (refs. 4 and 5) apply linear quadratic regu- 
lator theory to the design of controls for the F401 turbofan engine. 

The nonlinear system equations are linearized about five different 
equilibrium points , and linear system descriptions are obtained. The 
resulting linear models have five state and five control variables. 

At each equilibrium point, a quadratic performance index intended to 
minimize acceleration time is formulated, and feedback control gains 
are determined. A nonlinear feedback control law is developed by 
curvefitting the resulting control gains as a function of compressor 
speed. 

Weinberg (ref. 6) applies linear-quadratic regulator theory to 
the design of controls for the FIDO engine. He shows that this engine 
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can be adequately represented by three state variables - fan speed, 
compressor speed, and augmentor pressure. Four control variables are 
utxlized, and linearized engine models are obtained at two equilibrium 
points. The problem of minimizing acceleration time is considered, and 
control system gains are derived by conducting small perturbation opti- 
mizations at each of two equilibrium points, using a quadratic perform- 
ance index. The control gains are switched at a fixed value of fan 

speed, rather than varied in a continuous manner as in references 4 
and 5 . 

In reference 7, Sevich and Beattie consider the minimization of 
acceleration time for a turbojet engine, using fuel flow and exhaust 
nozzle area as control variables. They use a quadratic performance 
index to approximate a minimum time solution, as in references 4 to 6. 
However, they use a nonlinear engine model, rather than a series of 
linear models. The result is an open-loop optimal trajectory. The 
controls are assumed to be piecewise constant, and the performance 
index is minimized by using a conjugate gradient search technique. 

DeHoff et al, (ref. 8) use linear-quadratic theory to design con- 
trols for the FlOO engine. The control gains are generated using 
linear models with five state variables and four control variables at 
several equilibrium points. Principal emphasis is on the regulator 
design. Although acceleration control is considered, there is no 
specific attempt at minimizing acceleration time. 

References 4 through 8 all make use of integral quadratic per- 
formance indices, in which both state and control deviations from some 
desired trajectory are penalized. The coefficients of the penalty 
terms are adjusted in an attempt to minimize acceleration time. How- 
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ever, none of these reports can claim their final histories produce 
truly minimum-time accelerations. 

Minimum-time trajectory optimization has been considered by many 
investigators for many years. Athans and Falb (ref. 9) give a good 
account of the literature dealing with time optimal systems in their 
chapter 7. However, nearly all of their discussion is concerned with 
problems having only a single control variable. Furthermore, the sys- 
terns are assumed to be linear and time invariant, and the control 
limits are not dependent on the state or time. It is shown that the 
optimal control for such problems is bang-bang, i.e. , the control al- 
ways operates at either the upper or the lower limit. 

Wolske (ref. 10) considered the problem of fuel-optimal control 
of a dynamic system which is nonlinear in the state and linear in the 
control. The controls are assumed to be bounded in magnitude, and the 
resulting optimal control is bang-bang. The problem is solved by 
linearizing about a nominal history, which is neither optimal nor fea- 
sible (i.e., it does not satisfy the terminal constraints). The opti- 
mality condition and terminal constraints are expressed as linear in- 
equalities, and linear programming techniques are used to improve the 
solution until a feasible optimum is attained. 

C. Contributions 

In this report, an algorithm is developed for solution of optimal 
control problems which are nonlinear in the state variables , and linear 
in the control variables. Specifically, the problem considered is to 
minimize a performance index subject to satisfaction of the system 
dynamic equations , a set of terminal constraints (the number of which 
may be less than or equal to the number of states) and path inequality 
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constraints. The performance index and system equations and path con- 
straints are all linear in the control variables. 

It is shown that the optimal control is bang-bang,, except for 
possible singular arcs which are not considered. The algorithm re- 
quires that a nominal bang-bang solution be found that satisfies the 
system dynamic equations and terminal constraints. Once such a feasi- 
ble solution has been found , influence functions are generated which 
determine if the necessary conditions for optimality have been satis- 
fied. If not, additional control switches are added. Nonlinear opti- 
mization (gradient search) techniques are then used to vary the con- 
trol switching times in order to Improve the solution. The nonlinear 
optimization technique described in reference 11 is used to generate 
the numerical results presented in this report. 

The algorithm presented herein converts an optimal control problem 
with path inequality constraints and terminal constraints into an uncon' 
strained parameter optimization problem. This is accomplished in two 
steps. First, the bang-bang nature of the optimal control is used to 
express the possible optimal trajectories in terms of the switching 
times between regions of different control strategy. At this point, 
the original problem has been converted into a parameter optimization 
problem with equality constraints (terminal constraints) . Then, the 
equality constraints are satisfied by using an equal number of the 
switching times as iteration variables. This procedure results in an 
unconstrained parameter optimization problem with a reduced number of 
parameters, and is similar to the reduced gradient algorithm discussed 


in reference 12. 
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The algorithm is then, used to find minimum time acceleration his- 
tories for the FlOO engine, a two-spool turhofan engine used to power 
the F15 and F16 aircraft, A piecewise-linear engine model having three 
states and four controls is used to obtain the minimum time solutions. 
The linear models used in this report were obtained by linearization of 
a nonlinear model at five equilibrium points, and were taken from ref- 
erence 13. The linear model which applies at a given time in the tra- 
jectory is determined by calculating a normalized "distance" from the 
current state to the state at each of the equilibrium points; the 
linear model associated with the closest equilibrium point is then 
used. Linear state/control constraints which correspond to speed, tem- 
perature, pressure, and mechanical control limits are considered. Min- 
imum time solutions are obtained, and the resulting control histories 
are used as inputs to a nonlinear computer simulation of the FlOO 
engine (ref. 14) to verify the accuracy of the piecewise linear solu- 
tion. 

A suboptimal closed-loop control mode is also developed, which 
gives performance which closely approximates the open loop results. 

Use of the piecewise linear model allows optimal solutions to 
be obtained with the expenditure of less than one percent of the com- 
puter time which is required when a detailed, nonlinear model is used, 
such as in reference 7, Furthermore, the solutions obtained in this 
report are truly minimum time solutions to the piecewise linear prob- 
lem, rather than approximations to minimum time solutions to the non- 
linear problem, as in references 4 to 7, 

In addition to the transient results , the linear models and con- 


7 


straint equations are used to find the control settings which minimize 
steady-state specific fuel consumption « 

D. Organization of the Report 

The dynamic optimization problem is defined in chapter II* and 
necessary conditions for an optimal solution are stated. The algo- 
rithm is presented in chapter III. First, the initial feasible solu- 
tion is defined and discussed. Next, it is shown how to calculate the 
sensitivity functions (Lagrange multiplier functions) corresponding to 
the feasible solution. Finally, improvement of the feasible solution 
is discussed. Chapter IV derives the applicable equations for piece- 
wise linear systems. In chapter V, a detailed comparison of exact and 
piecewise linear solutions to a particular nonlinear problem is made. 
One and two control variable problems are discussed. Chapter VI pre- 
sents the results obtained for the minimization of acceleration time 
for the FlOO engine. Finally, conclusions are drawn and recommenda- 
tions for future research are made in chapter VII. 


I 
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CHAPTER II. A NONLINEAR OPTIMAL CONTROL PROBLEM, LINEAR IN CONTROL 
In this chapter, the optimal control of a dynamic system which is 
nonlinear in the state and linear in the control is considered. First, 
the problem is defined, including the performance criterion, system 
equations, path constraints, and terminal constraints. Then, necessary 
conditions for optimality are derived. The derivation of the optimal 
control strategy is similar to that presented in textbooks on modern 
optimal control, such as reference 15. This derivation is included 
here because the nature of the resulting optimal control strategy moti- 
vates the development of the new algorithm presented in chapter III. 

A. Problem Statement 

We consider a fairly general dynamic optimization problem, which 
is subject to one important restriction - the performance index, system 
equations and path constraints are all linear functions of the control 
variables. As will be shown later, this leads to a bang-bang solution. 
We wisix to find the vector control history u(t) which minimizes the 
scalar functional 


J = (j)[x(tp,t^] + 


[a(x,t) + b (x,t)u]dt (2.1) 


subject to the vector system differential equations 


X = f(x,t) + g(x,t)u 


( 2 . 2 ) 


and path inequality constraints 


c^(x,t) + d^(x,t)u <. 0 


X — 1,2, , . , , q 


(2.3) 
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The initial state and time are assumed to be specified, x.e. , 

x(tQ) - Xq 

while the terminal state and time are subject to the p terminal con- 
straints (p ^ n + 1) 

4 'i[x(tp,t£] =0 i = 1,2, . . .,p i (n + 1) (2.4) 

In the above, x is the (n x i) state vector, and u is the (r x 1) 
control vector. The functions <j), a, and c^ are scalar functxons of 
X and t, while b and d^ are vector functions of dimension (r x 1) . 
The vector function f and matrix function g have dimension (n x 1) 
and (n X r) , respectively. The terminal time t^ may be either fxxed 
or free. In fact, if a = b = 0 and = t^, the performance index in 

(2.1) is simply J = t^. 

The path constraints (2.3) serve to bound the allowable values of 

T 

the control. A path constraint is said to be active if + d^u - 0; 
it is said to be inactive if c . + d^u < 0 . If, for a particular path 
constraint, is constant and d^ has only one nonzero, constant 

component, then that path constraint is simply a physical control limit. 
On the other hand, if all components of d^ are equal to zero, the 
control does not appear explicitly in the path constraint} such con- 
straints are called state variable inequality constraints. In the main 
body of this report, it is assumed that there are no state variable in- 
equality constraints. However, the theory and numberical algorithm are 
extended to include state variable inequality constraints in appendix A 
If, for given values of x and t, the component of d^(x,t) 

is positive, then the i^^ constraint serves as an upper bound for the 
control variable. Similarly, if the j component of dj^(x,t) is 
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negative, the constraint is a lower bound for the control var- 
iable, It will be assumed that sufficient constraints are imposed so 
that each of the r control variables is bounded above and below for 
all possible values of x and t. If this is the case, then nc con- 
trol impulses are allowed. 

B. Necessary Conditions for a Local Mnimum 
Using the techniques employed in reference 15, a Hamiltonian func- 
tion is defined as 

H ^ a(x,t) + b^(x,t)u + A^[f(x,t) + g(x,t)u] 

q 

(x,t)u] (2.5) 

where X and p are undetermined Lagrange multiplier vector functions 
of time, having dimension (n x 1) and (q x 1) , respectively. 

In terms of H, necessary conditions for J to be a local minimum 

are 


and 

The Lagrange multipliers must satisfy the following teirminal con- 
ditions. 

|i^[x(tp,t£l + ||- [xCtp.Cf] (2.7) 


•T _ 3H 
^ ” 3x 


3u 


= 0 


(2.6a) 

(2.6b) 


if c, + d.u = 0 
i X 


(2.6c) 


if c. + d.u < 0 

X 1 


^j^[c^(x,t) + d^ 
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where v is a (p x 1) undetermined parameter vector. If the terminal 
time is not specified, another necessary condition which applies is 

3(|)[x(t^) ,t^] ^ 94'[x(tj) ,tj] 


H(tp = - 


9t, 


- V 


9t, 


( 2 . 8 ) 


■f “"f 

For the present problem, equations (2.6a, b, and c) may be written 
explicitly as 

T ^ I 

9a 9b^ ,T 9f .T ^ i q. 


i=l 


T T 
b^ + X^g + 


D = 

i=l 


(2.9b) 


^i 

= 0 

if 

c , . + 

X 

d7u 

X 

< 0 

^i 

|v 

o 

if 

c. + 

X 

d^u 

= 0 


(2.9c) 


The optimal control must satisfy (2.9b and c) and constraints 
(2.3). For given x and t, the functions b, c, d, g, and X are 
all constant. Therefore, the determination of the optimal u (for 
each X and t) is simply a linear programming problem, which can be 
readily solved by using the Simplex method (ref? 16) . Except for 
singular arcs (to be discussed shortly) , the optimal control is always 
determined by r active constraints from (2.3). The optimization 
problem is simply to find the r constraints which satisfy (2.9b and c) 
and (2.3). 

The r active constraints change from time to time along the tra- 
jectory. When a change of constraints occurs, the control variables 
jump dis continuously from one boundary to another. Such control is re- 
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f erred to as bang-bang control, and the points at which the jumps 
occur are called junction points. 

C. Singular Arcs 

The optimal control can usually uniquely determined from (2,9b 
and c) and (2.3). However, if H is constant along an active con- 
straint boundary, i.e.,if b + Ag = where n is a real vari- 

able, then the control is not uniquely determined along that constraint 
boundary. This situation is illustrated in sketch (a) for a two con- 
trol variable problem. 



increasing H 


Iiir the sketch, H is constant along constraint and it is not clear 
whether to use constraint © or constraint ^^^(or neither) as the 
other active Constraint. Furthermore, if (b^ + A^g) = 0 (all compo- 
nents zero) the control is totally indeterminate. 

It is sometimes possible to -find minimizing solutions for which 
some or all of the controls u are not determinate from (2.9b and c) 
and (2.3) for a finite time period; the corresponding trajectory seg- 
ments are called singular arcs . On singular arcs, the control is de- 
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termined from the requirement that the time derivatxves of 
(b^ + X^g - ndj). as well as (b'^ + X^g - ndj) itself, must be identi- 
cally zero. To determine the control u, successive time derivatives 
of (b^ + X^g - ndT) are taken until the control appears explicitly. 

It is difficult to determine general conditions for the existence 


and minimality of singular arcs. Nevertheless, each individual problem 
should be examined for the possibility of minimizing singular arcs, and 
this practice will be followed in this report in the example problem to 
be discussed later. However, in the algorithm to be described in chap- 
ter III, it will be assumed that the minimizing solution is nonsingular, 
p. Determination of Optimal Trajectory 
We now assume that the problem is nonsingular, and the linear pro- 
gramming problem is solved to yield the active constraints and the 
optimal control. If we assume that the active constraints are con- 
straints 1 through r and form these into a vector-matrix represen- 
tation, i.e.. 


C = 




then we can solve for the optimal control from 


-T 

u = -D C 


( 2 . 10 ) 


where the superscript (-T) denotes the inverse of the transpose. Also, 
equation (2.9b) can be rewritten using vector-matrix representation 


«-Via constraints 
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CHAPTER III. A NEW ALGORITHM 

Necessary conditions for the optimal control of a dynamic system 
which is nonlinear in the state and linear in the control were derived 
in Chapter II. It was shown that except for singular arcs, the opti- 
mal control (r variables) is always determined by r active con- 
straint boundaries. In Chapter III, the nature of the optimal control 
strategy derived in Chapter II is used as the basis for a new" algorithm 
for the solution of such optimization problems. 

a feasible solution is obtained which satisfies all path 
and terminal constraints , and for which the controls always lie along 
r constraint boundaries. The Euler Lagrange equations are not uti- 
in the determinatxon of this feasible solutxonj it may or may 
not be a local minimum. Then, it is shown that the Lagrange multiplier 
time history can be easily and uniquely calculated from this feasible 
solution. The necessary conditions for a local minimum may be calcu- 
lated as functions of the Lagrange multipliers. If the initial fea- 
sible solution does not satisfy the necessary conditions, the control 
history is modified, and nonlinear optimization (gradient search) tech- 
niques are used to improve the solution. 

A, Modified Problem Statement 

In chapter II it was shown that except for singular arcs, the 
optimal control variables (r variables) are always determined by an 
equal number of active constraints at all points along an optimal tra- 
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jectory. The problem is to find which constraints are active, as a 
function of tine. It will be assumed that the optimal trajectory con~ 
sists of a number of segments. For each segment, the control variables 
are determined by a set of r active constraints. The optimal control 
problem stated in (2.1) to (2.4) is solved by finding the optimal 
values of the switching times, i.e. , the times at which the active 
constraints switch from one set to another. 

The modified optimal control problem can be stated as follows: It 

is desired to find the values of the w switching times , 

{ t^ , i = 1 , . . . , w } 

which minimize the scalar functional 


J = cf)[x(tp,tj] + 



[a(x,t) + b (x,t)u]dt 


while satisfying the vector system differential equations 

X = f(x,t) + g(x,t)u 

and terminal constraints 

ii^^[x(tp,t^] =0 i - 1, . . ., p < (n + 1) 


During time interval 


the control is given by 


t. < t ^ t 
X-1 X 


u 



C. 

X 


in accordance with (2.10). 

The above modified problem is a parameter optimization problem 
with equality constraints, which might be easier to solve than the 
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original optimal control problem. However, it should be emphasized 
that the number of segments in an optimal trajectory, and the active 
constraint sets for each of the segments, are not known in advance. 

The problem statement can be further modified into an uncon- 
strained parameter optimization problem by using the reduced gradient 
method, as presented in reference 12. To accomplish this, the switch- 
ing time set is partitioned into two subsets, having p and (w - p) 
members, respectively. The p members of the first subset are con- 
'sidered as dependent variables, in the following sense: Whenever some 

of the switching times in the second subset are varied, the values of 
the p terminal constraints will change from their converged values. 
The values of the p switching times in the first subset are then 
used as iteration variables in order to reconverge the terminal con- 
straints to the desired final values. In this way, the problem is 
converted into the following form; Find the values of the (w - p) 
switching times 

i = (p + 1), . . • , w) 

which minimize the scalar functional 
J = ^[x(tp,tj] + 

while satisfying the vector system differential equations 

X = f(x,t) + g(x,t)u 



[a(x,t) + bHx,t)u]dt 


where, during time interval 


i I f 
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the control is given by 


-T 

u = -D. C. 

X 1 

The above unconstrained parameter optimization problem may be solved 
by using any of a number of nonlinear optimization techniques avail- 
able, such as those discussed in reference 12. The results presented 
in this report were obtained by using the parameter optimization tech- 
nique reported in reference 11, 

B. Feasible Solution 

In order for a trajectory to be a local minimum solution to (2.1), 
equations (2.12) and (2,13) of Chapter II must be satisfied, and the 
control must satisfy constraints (2.3) and optimality conditions (2.9b 
and c). In addition, terminal constraints (2.4) and (2.7) must be 
satisfied, and equation (2.8) must be satisfied if the terminal time 
tf is free. In developing the algorithm for the solution of this 
problem, it will be assumed initially that t^ is free. Later, we 

will show how to modify the algorithm for the case in which t^ is 
specified. 

We define a feasible solution as one which satisfies the system 
differential equations (2.2) and the terminal constraints (2.4), and 
where the control u(t) is consistent with the necessary conditions 
for optimality - that is, the control is determined by r of the 
path constraints (2.3) at all points along the trajectory. It is not 
necessary that the control be determined by the same r constraints 
at all points along the trajectory - in fact, we will usually require 
the control to be determined by several different constraint sets, as 
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will be seen shortly. Such a trajectory may or may not be a local min- 
imum solution for the performance index (2.1) o 

In order to obtain a feasible solution, there are p terminal 
constraints (eqsr (2.4)) which must be satisfied. In general, there 
must be p degrees of freedom available in order to satisfy the p 
terminal constraints. In order to provide these degrees of freedom, 
it will be assumed that there are at least p segments in the trajec- 
tory. For each segment, the control is determined by choosing the r 
constraints which are active. The set of active constraints may be 
chosen arbitrarily for the first segment; for each succeeding segment, 
one of the active constraints should be different than any which was 
utilized in the preceding segments. The durations of p of the seg- 
ments are variable, and provide the p degrees of freedom necessary 
to satisfy the p terminal constraints. 

The choice of the constraints to be active for the various tra- 
jectory segments should be made carefully. There is no guarantee that 
a solution exists for arbitrary choice of the active constraints, or 
for any choice of active constraints, for that matter. Also, some of 
the constraints assumed to be inactive in the initial feasible solution 
may be violated. In this case, different active constraints must be 
chosen. 

The combination of p degrees of freedom and p terminal con- 
straints is known as a multipoint boundary value problem, and must gen- 
erally be solved iteratively. Two widely used classes of methods for 
solution of such problems are Newton-Raphson methods and gradient 
methods. Both of these methods are widely discussed in the literature 


20 


(see, for example, refs. 12 and 17) and will not be discussed here. 
However, both methods require initial guesses for the values of the p 
degrees of freedom, and partial derivatives of the terminal constraints 
with respect to the degrees of freedom. Equations for calculation of 
the required partial derivatives are derived below. 

Suppose, for example, that initial guesses are made for (p - 1) 
of the junction times t^, and for t^. The choice of these times is 
sufficient to determine a reference trajectory, and values of the p 


terminal constraints ilj,. 

d 

If one of the t. is altered slightly, the state at t, is 


changed by 


dx(t.) = [x(t ) - x(t1‘)]dt. 


The effect of changes in the t, on the terminal constraints 


may 


therefore be obtained by integration of 


A. (al\ - 

dt \hJ - 


4 If: 


i = 1, 2, . . . , (p - 1) 


from t - t^ to t - t^, with initial conditions 


(3.1a) 


1 


(3.1b) 


The 9i]j^/9t^ are then calculated from 


. 9x 


3t . 9x_ 9t . 

X f X 


(3.2a) 


9i|;. 9i|j, 

9t^ 9Xj 


x(t.) 


(3,2b) 


Also, we have 


fJk A 
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I jjced final time . - it was assumed in the above discussion that 
the terminal time t^ is free. If t^ is fixed, then the duration 
of the final segment does not provide a degree of freedom. Therefore, 
there must be at least (p + 1) segments for p terminal constraints, 

and the durations of p of the segments provide the necessary degrees 
of freedom. 


C. Calculation of Lagrange Multipliers 
We consider first the case in which the terminal time is free. 
Suppose a feasible solution has been found. We will show that the 
Lagrange multipliers A and y can be uniquely calculated, as a 
function of time. Once the multipliers have been calculated, the 
necessary conditions for optimality can be checked. If the necessary 
conditions are not satisfied, an iterative improvement scheme is used 
in order to find a local minimum solution. 

Let the i^^ trajectory segment have initial time t^_j^ and final 
time t^. The control for the i*^^ segment was determined by a set of 
r active constraints, and the control just prior to t^ is denoted 
by u(t^). Just after t = t., the control is determined by different 
active constraints, and is given by u(t^) ^ u(t~) . it is shown in 

appendix B that the Hamiltonian must be continuous at t = t^. There- 
fore, we must have 


b 


b (t . ) + A^g(t 


.^|u(t+) - u(t-^ =0 1 = 1, 2, . . 


(P - 1) 


(3.3) 


where the notation g(t^) is shorthand for 

S(t^) ^ g[x(t^) ,t ] 
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Equation (3.3) must be satisfied at each of the (p - 1) junction 
points which are varied in order to satisfy the terminal constraints. 
Also, as shown in appendix B, the terminal X and H must be given 
by (2.7) and (2.8). 


T,, . T H 


^^^f^ " 3t 


-V 


at, 


If we use (2.5) for H, and substitute for u and A from (2.10) and 
(2.7), respectively, the above equation is converted to 


T dijj 
^ dt 


(t.) 


dt^(t^) 

+ - + a(tp + b^(tpD (tpc(tj) = 0 


(3.4) 


Equations (3.3) and (3.4) give p equations which must be satisfied, 
and there are p multipliers which may be varied. Thus, we have a 


multipoint boundary value problem, similar to that which must be solved 
iteratively to determine a feasible trajectory. However, in the pres- 
ent case, iterative solution is not required. Instead we can solve 
for the parameters v as follows: 

First, we find (p + 1) backward solutions of the A equation 
(2.13) for = (0, . . .,0), v'^^(l= 0, . . .,0), 

= (0, 1, 0, . . ., 0), . . ., v'^= (0, 0, . . ., 0, 1) where 


X^(t^) 


= 4- J 

8x^ ^ ^ ax^ 


These backward solutions are called 

A*^°\t), A^^^(t), . . . , A^P^(t) 


(3.5) 
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Then A.(t) is given by 

X(t) = + A(t)v 

where A(t) is defined as 

A(t) k (t) • . . . • X^(t)] 

• • 

There are (p - 1) equations for the continuity of H: 


[u(t^) - u(tp]^[b(t.) + g^(t.)X^°^(t.) + g"^(t.)A(t.)v] = 0 

i = 1, 2, . . . , (p - 1) 

and at t^, we must have 


a(tp - b^(tpD ^(tpc(tp (tp + 


d(|) 


f ' 


V = 0 


These equations may be put in matrix form as follows; 

Define vectors r. and q. and matrices Q and R by 


1 1 

T A r ,.+v ..-mT _T 


qT k [u(tp - g^(t^)A(t^) i = 1, .... (p - D 

rj i [u(tt) - u(t")]’^[b(tp + g'^(t^)A^°\(tpl 

T i f 


s - 


dt 


:l k a(tj) +b'^(tpD-’^(tpo(tp +f| (tj) 


Q 4 



/ 


R k 


Then V may be calculated from 


(3.6) 

(3.7) 

(3.8) 

(3.9) 

(3.10) 

(3.11) 


1 


1 
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V = -Q“^R (3.12) 

T 

Once the v are know" ^ (2.7) can be used to calculate A (t^) , 
and A(t) can be obtained by integrating (2.13) backward in time, or 
by using (3.6). 

Fixed final time . - In the above development, it was assumed that 
the terminal time t^ is free. In the event that t^ is fixed, the 
equations for calculation of v are easily modified. In this case, 
(3.4) is not applicable. Instead, there are p of equations (3.8), 
instead of (p - 1) , since there are (p + 1) segments for this case. 
The p equations (3.8) are sufficient to calculate the p param- 
eters, V. 

D. Improvement of Feasible Solution 
The Lagrange multipliers can be used to determine if the neces- 
sary conditions for optimality are satisfied by the initial feasible 
solution. The optimal control is obtained as a function of time by 
finding active constraints which result in satisfying (2.9b and c) and 
(2.3), using x(t) and A(t) as determined from the feasible solution. 

If u . (t) as determined in this manner is identical to the control 
opt 

time history u^^(t) utilized in the feasible solution, then the ini- 
tial feasible solution satisfies the necessary conditions for an opti- 
mal solution, and no further calculations need be made. On the other 
hand, if u^^^(t) differs from u^^(t) for even a portion of the tra- 
jectory, then the feasible solution is not a local minimum. 


Suppose, for example, the control 


u . (t) differs 
opt 


from 




during trajectory segment k. Then the performance index (2.1) can be 
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improved by splitting segment k into two parts , and using control 
as follows . 


u(t) = Uj^(t) , i t 1 


(3.13) 


where t^^ should be chosen to be only slightly greater than t^ 
so that the modified trajectory differs only slightly from the initial 
feasible trajectory. Because of’ the modified control history, the new 
trajectory will not satisfy the p terminal constraints (2.4). There- 
fore, the original p junction times should be adjusted, while hold- 
ing t fixed, so that the terminal constraints are satisfied. 

During the time interval from to t^^, the control has 

been changed (from that of the initial feasible trajectory) from 
u^^(t) to u^p^(t). Therefore, the state at t^^ is changed (to 
first order) by 

Ax = g(t ) [u (t ) - u- (t ) ] (t - t, - ) (3.14) 

sw * s\j opt sw^ fe sw sw k-1 

The effect of this change on the terminal conditions is deter^ 

mined from 




Ax(tp 


(3.15) 


where Ax(tp is obtained by integration of 

~ (Ax) = (f - gD“^C)(Ax) (3.16) 

from t = t^^ to t - t^, with initial conditions given by (3.14). 

The change in the must be canceled by appropriate changes 

in the t.. Therefore, we must have 
1 ’ 
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9ii;. 

AiJj . + At . = 0 
1 9t. X 

X 

which may be solved for At^ to give 



Equation (3.17) is applied repeatedly in an iterative manner until the 
terminal constraints are satisfied to sufficient accuracy so that fur- 
ther iteration is not required. 

Once a modified trajectory has been obtained and the -terminal con- 
straints satisfied, the Lagrange multiplier time history may be calcu- 
lated for the modified trajectory in the same manner as it was calcu- 
lated for the initial feasible trajectory. From the Lagrange multi- 
pliers, the gradient of the performance index with respect to the 

switching time, i.e. , 9J/3t can be obtained. It is shown in 

sw 

appendix B that 

H(t" ) - HCt"^ ) = A^(t )g(t )[u(t" ) - nit*' )] 

9t sw sw sw sw sw"^ sw 

sw 

+ b’^(t ) [u(t“ ) - uCt"^ )] (3.18) 

sw sw sw 

The set of all such switching points t^^ and corresponding 
gradients, can be used in conjunction with a nonlinear search tech- 
nique (see refs. 12 or 17 for a general discussion of nonlinear 
search techniques, or ref. .11 for the particular search technique used 
to obtain the results presented in this report) to search for the 
values of t such that the gradient vector 3J/3t is equal to, 
or nearly equal to, zero. 
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E. Merging of Switch Points 

In this report, optimal solutions are found by searching for the 
values of a number of control switching times which minimize a perform- 
ance index and satisfy terminal constraints. These switching times 
divide regions of different control action, in which different path in- 
equality constraints are active. In addition to the number of such 
control regions, the order in which these regions occur (timewise) must 
be specified initially. However, during the course of the search for 
an optimal solution, the switching times do not necessarily remain in 
the originally specified order. 

Suppose, for example, the set of switching times for the initial 

feasible solution is given by {t^, i = 1 w ^ p> and 

< tj^ < . . . < t^ < tj. The values of p of the switching times 
are varied in order to satisfy the terminal constraints, while the re- 
maining (w - p) switching times are varied to minimize the performance 
index. During the search for an optimal solution, there are three 
possible situations wliich can arise which affect the order and/or num- 
ber of switch points. These are: (1) t^ becomes less than t^ ; 

(2) t^ becomes greater than t^; (3) t^ becomes greater than 
where 1 i k £ (w - 1) , 

Cases (1) and (2) may be handled in an identical manner. It is 
tentatively assumed that segment 1 (or segment (w +1)) is not required, 
and the search is restarted with one less segment. Once this reduced- 
segment search has converged, the necessary conditions for optimality 
are checked, and the segment which had been eliminated is reinstituted 
if necessary (along with other additional segments, as appropriate). 


The segment which was eliminated may have been one of the p segments 
used to satisfy the terminal constraints. In this case, one of the 
remaining segments which had been used to minimize the performance 
index must be used instead to satisfy the terminal constraints. 

The procedure followed if case (3) is encountered involves several 
possibilities, and can result in a change of active constraints or a 
loss of one or two switching times. It was stated earlier that at each 
switching time, one of the active constraints changes. For definite- 


ness, let the active constraint from segment k which is inactive in 
segment (k + 1) be denoted by constraint a. Also, let constraint b 
denote the new active constraint for segment (k 4- 1) . 

Similarly, let constraint c denote the active constraint from 
segment (k + 1) which is inactive in segment (k + 2) , and let con- 
straint d denote the new active constraint in segment (k + 2) . 
Clearly, a?^b, c?^d, a?^c, and b =/ d, but b and c, and/or a 
and d, may be the same. The several possibilities are as follows: 

(3.1) b c, a d. For this case, if t, and t, become 

k k+1 

equal and interchange, then the active constraints for segment (k + 1) 
become a and d, instead of b and c. No segments are eliminated. 

(3.2) b = c, a d. For this case, if t, and t, become 

k k+1 

equal, segment (k + 1) and switching time eliminated. 

(3.3) b = c, a = d. For this case, if t. and t, become 

k k+1 

equal, segments k and (k + 2) merge into one, and segments k and 
(k + 1) , and Switching times tj^ and t^^^ are both eliminated. 
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F. Insertion of New Switching Times 
Suppose a local minimum solution has been found for a particular 
set of switching times {t^, i = 1, » . w ^p}« As discussed earlier, 
the next step is to determine if the necessary conditions for a local 
minimum solution have been satisfied. If not, new switching times are 
added, and a new search for a local minimum is undertaken. The new 
switching times are added in the region in which the necessary condi- 
tions are violated. Initially, the new switching times are inserted 
at times exactly equal to existing switching times. In this way, new 
segments having new active constraints are added, but with zero time 
duration. 

For example, suppose an initial feasible solution is found having 
w segments , and it is discovered that the necessary conditions for a 
local minimum are not satisfied during segment k. Then a new segment 
and corresponding switching time are added during segment k. Initi- 
ally, the new switching time is placed at = t^_^ so that the 

new segment starts at tj^ ^ and has zero duration. A segment having 
zero duration cannot alter the trajectory or performance index; there- 
fore, we have 

j(t^, . . ., t^, = Vl^ " ^ *^w^ 

i.e., the performance index is continuous with respect to insertion 
of new switching times. 

G. Convergence Issues 

It is highly desirable that one be able to demonstrate in ad- 
vance that a particular search will always converge to a solution 
which satisfies the necessary conditions for a local minimum. There 
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are several properties which the algorithm must have in order to 
assure such convergence. First, the algorithm must he setup as a 
descent algorithm. That is, each trial solution must have a smaller 
performance index than the preceding solution. This descent property 
is possessed by most available nonlinear optimization techniques, in- 
cluding the one used herein (ref. 11). Additionally, both the algo- 
rithm and the performance index riust be continuous functions of thexr 

arguments, in this case the switching times. 

It is clear that the performance index is a continuous function 
of the switching times so long as the switching times do not inter 
sect, and no switching points are deleted or inserted. It remains to 
be shown that the performance index is continuous even when switching 
times intersect, and/or switching times are added or deleted. 

When switching times intersect, the algorithmic search is termi- 
nated. As discussed previously, the intersection may result in the 
deletion of zero, one or two switching points. In any case, the de- 
leted switching times correspond to segments having zero duration at 
the time of intersection. Since segments of zero duration do not 
affect the trajectory or performance index, it can be concluded that 
the performance index is continuous with respect to switching time 
deletions. After the deletion of switching times, the search is re- 
started with the appropriate number of switching times. 

Switching times may also be added after a search is completed if 
it is discovered that the necessary conditions for a local minimum 
solution have not been satisfied. This possibility has been discussed 
previously. It was shown that the performance index is continuous 
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i 


with respect to such switching time additions because the new switch- 
ing txmes are added in such a manner as to create segments of zero 
duration initially. 

The one required property which unfortunately is not possessed 
is continuity of the algorithm. Whenever switching times are added 
or deleted, the number of arguments changes, and it can be shown that 
the algorithm is discontinuous under such conditions. Because of this, 
convergence to a solution which satisfies the necessary conditions for 
a local minimum cannot be guaranteedc Nevertheless, such convergence 
is assured if no switching time insertions or deletions are encoun- 
tered. Furthermore, convergence to a solution which satisfies the 
necessary conditions for a local minimum will be attained in most in- 
stances even if switching time insertions or deletions are encountered 
during the search. 

H. Summary of Algorithm Steps 

It is useful to summarize the steps which are followed in the 
determination of an optimal solution, using the algorithm presented 
herein; 

(1) An initial feasible solution must be found. The number of 
segments w^p is selected, and the r active constraints for each 
segment are chosen. The p segments which are varied in order to 
satisfy the p terminal constiatats are alsc? chosen. The initial 
feasible solution is obtained by iterating on the p variable switch- 
ing times until the p terminal constraints are satisfied to desired 
accuracy. (2) The Lagrange multipliers corresponding to the initial 
feasible solution are calculated by using the procedure described in 
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Chapter III. The necessary conditions for optimality are calculated 
from functions of the Lagrange multipliers. (3) If the necessary con- 
ditions for optimality are not satisfied, additional control segments 
are added (if required) and the values of the („ - p) switching times 
are used in conjunction with a parameter optimisation (gradient search) 
procedure to improve the initial feasible solution. 


CHAPTER IV. PIECEWISE-LINEAR MODELS 


An important special case of the problem considered in chap- 
ters II and III is when the performance index, system equations and 
path constraints are all linear in both the state and the control. 

This case occurs frequently in practice because linear approximations 
to complex dynamic systems are readily available. Furthermore, solu- 
tions to linear problems are easier and less costly to obtain because 
the system and Euler-La grange equations may be represented by transi- 
tion matrices , 

If the actual system is only slightly nonlinear, a single linear 
approximation to the nonlinear system may suffice over the full oper- 
ating range. However, if greater accuracy is desired and/or the 
actual system is very nonlinear, a series of linear models may be used, 
each of which is obtained by linearizing about a different equilibrium 
point. Linear equations are still used to describe the system at each 
state point, but the coefficients in the linear model vary from point 
to point. Such a model is called a piecewise linear model. 

In this chapter, linear and piecewise linear models will be con- 
sidered. Although the development of chapters II and III is fully 
applicable to this problem, significant results for the linear problem 
will be repeated here because of the importance of this special case. 
The one-piece linear model will be considered first; then, the piece- 
wise linear model will be considered. 
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A. Linear Model 

The problem to be considered is to find the control u(t) which 
minimizes the functional 

nh 

J = (}i[x(tp,t^] + / (a^x + b"^u)dt (4.1) 

subject to the system differential equations 


X = Tx + Gu + h 


(4.2) 


and path inequality constraints 


cTx + dTu + e.iO i = l,2,. . .,q 

X 1 X ’ ’ ’ ^ 


(4.3) 


Here, the vectors a, b, c^, d^, and h, and the matrices T and G, 
are all constant. The terminal constraints on the state are given by 
(2.4). The Hamiltonian for the linear problem is 

q 


T 

H - a X 


+ b"^u + l^(Fx + Gu + h) + ^^V^(cTx + d^u.+ e^X^« 

i=l 


4) 


and the resulting Euler-Lagrange equations are 

q 


A. = -a - A F 


- 

i=l 


(4.5) 


The Lagrange, multipliers must satisfy terminal conditions (2.7) and 
(2.8). The optimal control is determined from 

q 


- 0 

i=l 


(4.6) 


where 


35 


> 0 

if 

i 

c,x 

X 

+ 

d'u 

X 

+ e . 
X 

= 0 

= 0 

if 

T 

c,x 

1 

+ 

d .u 

X 

+ e, 

X 

< 0 


We assume, as in chapter II, that the solution is nonsingular 
and the active constraints and optimal control have been determined 
from (4.6) and (2.3). The active constraints are formed into matrices 
as follows : 


C = 


/ T \ 

\ 

\ T / 

\ r / 



and the optimal control is given by 


-T 

u = -D Cx 


The multipliers y given by 


y = -D“^(b + G^A) 


(4.7) 


(4.8) 


Substitution of (4.7) and (4.8) into (4.2) and (4.5) results in 


X = (F - GD ^C)x + h (4.9) 

A = -(a^ - b^D"-^C) - A (F - GD C) (4.10) 

Determination of the feasible trajectory and Lagrange multipliers 
proceeds exactly as in chapter III. However, the sensitivity 
vectors defined in (3.1) and (3.5) as well as solutions to the system 
equations and Euler-Lagrange equations, may be represented by transi- 
tion matrices because the equations are linear. For example, (3.1) 
becomes 


J 


I 


I 


dt 


(%)■ 
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(F - GD"^C) 


te) 


which has the general form 

^=M+NP; F(0) = Pq (4.11) 

where P and M are (n x 1) vectors and N is an (n x n) matrix. 

Tlie solution of (4.11) is given by 


P(t) = <S>(t)PQ + 



$(t - t)M dt 


(4.12) 


or 


P(t) = ^(t)PQ + N"^[$(t) - I]M 


where 

(J> = 

with initial condition 

$(0) == I 

Also, $(t) may be calculated from 


$(t) = TA(t)T' 


where 


A(t) = diag 


A- is the i*"^ eigenvalue of N, and the i^*^ column of T is the i*"^ 
eigenvector of N. 

B. Piecewise-Linear Modeling 

The desirability of using linear equations to model a system is 
obvious. Nevertheless, the actual system may be so nonlinear that a 
linear system description is not sufficiently accurate over the full 


A 


^nt] 

. . . , e J 

_.th , 
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range of interest. In such cases, it is possible to increase model 
accuracy while retaining the advantages of linear modeling by using a 
piecewise-linear system model. 

Suppose, for example, the nonlinear system is linearized about a 
number of equilibrium points. In the neighborhood of each equilibrium 


point, we have 


X = F^(x - x^^) + G^(u - u^^) 


( 4 . 13 ) 


where x . and u^^ are the equilibrium values of state and control 
at the equilibrium point j. The system matrices and also 

differ, in general, for each equilibrium point. The path constraints 
may also be linearized about each equilibrium point to yield 


'*11 + ®lj 




S3 


The path constraint vectors c^^ and d_ also differ for each equi- 
librium point. 

With a piecewise linear model, the system is described by linear 
equations at each state point, but the linear system coefficients vary 
from one point to another. A question that arises is which equilibrium 
model applies best for a given state, x? It is natural to choose the 
equilibrium point which is closest to x in some sense. Since the 
various states do not necessarily have the same physical dimension, a 
normalized distance function is used to determine which equilibrium 
point is closest. For a given state x, the distance functions 


I. = (X - x^.) W(x - x^.) 


( 4 . 15 ) 


are calculated for each equilibrium point j, and the equilibrium 
point is chosen for which I. is a minimum. 
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Since the I. are continuous functions of x, model switches 
3 

(say, from j to k) occur when I. = I, or 

3 ** 


(x - - x^^) = (x - X^^,)\(x - X^J^) 

Simplification of this expression reinilts in 


If we define 


o/ T T T T 

2(x 1 - X ,)Wx - X , Wx , - X ,Wx , 

^ ek ej ek ek ej ej 




where 


V - 


(4.16) 


(4.17) 


and 


ejk - 4«^ek - 

then, switches between equilibrium models j and k occur at 

S.i ^ 0 

3^ 


(4.18) 


Although (4.16) is linear in x, iterative solution is generally re- 
quired to find the switching points , since x is not a linear func- 
tion of time. 

C. Necessary Conditions for Optimality 
The problem to be solved is to minimize 

J s= ^[x(t^),t^] + / 

^^0 


(a^^x + b"^u)dt 


subject to the system differential equations 


x=F.(x-x.) +G.(u-u .) 
3 2 &3 
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and path inequality constraints 


Cij(x - + d^.Cu - u^.) + iO 

i— 1,2, . • »,q 

where the linear model index j is determined from 

T 

min(x - X ) W(x - x ) 

3 ® J *'3 

At switching points between models j and k, we have 

T 

3 ^ 3k 3k 

where a,, and B., are defined by (4.17). The Hamiltonian is 
jk jk. 

defined as 

H = a"^x + b^u + x'^LF-Cx - x ,) + G . (u - u .)] 

3 3 *=3 


i=l 


(X - x^j) + d.^(u - u^.) + e^.l 


and the Euler-Lag range equations are 




4 

VA T 

j - 

i*! 


where the model index 3 is the same function of time as determined 
by integration of the system equations. In addition, it is shown in 
reference 15 that the Lagrange multipliers are discontinuous at model 
switching points, the jump in X being given by 


X(t^) = X(t;) + 


(4.19) 


where t is the time at which the model subscript changes from j 
s 

to k. Furthermore, since the switching time is not specified a priori, 
the Hamiltonian must be continuous at t = t^. Therefore, we must have 


AO 


b’'u(t+) + [>.(0 + 


Jk 


(4.20) 


which gives 


e = 


x'^(t;)[x(t;) - x(4)] +b'^[u(t;) - u(e^)i 


T ‘f^+\ 


(4.21) 


The equations (4.6) for determination of u, the equations (4.7) 
and (4.8) for u and u, and the equations (4.9) and (4.10) for the 
working forms of the system and Euler-Lagrange equations, all apply 
to the piecewise-linear model, but with the appropriate model sub- 
script 3 added. 

D. Initial Feasible Trajectory 

Calculation of the Initial feasible trajectory proceeds in the 

saine manner as in chapter III, except that the model switches must be 

made at the appropriate times. Also, the calculation of the partial 

derivatives 3it.,/3tj^ in (3.2) must be modified. Suppose the Junction 

time t. is perturbed slightly and there is a model switching time 

t > t . Since t is determined by satisfying (4.18), the change in 
’"s i* s 

t. causes t^ to change. Therefore, the change in t. not only 
aHects the * directly as in (3.2), but also indirectly through the 

ef fect of a change in on • 

Specifically , a change in causes x(tg) to change by 

3x_ 


1 


( 4 . 22 ) 


Where 3Xg/3t. is obtained by integration of (3.1) from t. to t^ 
The resulting change in the model switching functxon S is 
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T T 

dS = a dx(t ) = a"^ -r— - dt. 

S dt . X 


This change in S must be canceled by a change in t^. 


we have 


T T. — 

dS = a r— dt. + a x(t )dt = 0 
ot. X s s 

X 


from which it follows that 


T 


at 

s 

8t. 

X 


at. 

X 


a^x(t“) 

s' 


Therefore, ai|)./at^ is given by 


3 

3ii;. 
= 3 


axf 

at. 

3x 

at. 

at 

X 

f 

X 

s 


T ^^s 

* at. 

X 


(4.23) 

Therefore, 


(4.24) 


(4.25) 


E. Calculation of Lagrange Multipliers 
The procedure for calculating the Lagrange multipliers correspond- 
ing to the initial feasible solution differs from that employed in 


chapter III, because of the jumps which occur in the multipliers at 
model switching times, as given by (4.19) to (4.21). For simplicity, 
we consider the case in which there is a single model switching point. 
The procedure which will be derived can be extended to the case of mul- 
tiple model switches. At a model switching point, we have 

^’'(0 - (4.19) 

As in chapter III, the value of H must be continuous at the (p - 1) 
variable junction points, t = t^, . , ., t^_^. In addition, for final 
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time free, equation (3.4) must also be satisfied. 

Except at t , the requirement for continuity of H is achieved 
s 

by satisfying (3.3). The requirement for continuity of H at t = 
requires special consideration because of the change in the system equa- 
tions , and can be expressed as 


AH = b^[u(tg) - u(tpi + X^(tg)[i(tg) - i(tg)] + ea^^x(tg) (4.26) 

In order to solve for e and the p values of v, we find 
(p + 2) backward soluciona of the A equation (2.13). The first 
(p + 1) of these solutions are identical to (3.5); (t) is ob- 

tained by integrating (2.13) backward, starting at t = t^, with ini- 
tial conditions 




By superposition, A (t) is given by 

A(t) = A^°^(t) + A(t)v + eA^P'*'^\t) (4.27) 

The (p + 1) equations for the determination of e and v are given by 

[u(t'l') - u(t7)]^[b(t ) + g'^(t )A^°\t ) + g^(t )A(t )v 

X X X -L 

+ g^(t.)A^P-^^>(t)e] =0 i = 1, . . .. (P - 1) (4.28) 

[u(tg) - u(tp]'^b(tg) + [x(tg) - x(t“)]'^[A^°^(tg) + A(tg)v] 

+ x(t'^)^a.,e = 0 (4.29) 

s 3 ic 


a(t^) - b^(tpD‘”^(tpc(tp 




(4.30) 
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The calculation of v and e proceeds exactly as in (3.10) to (3.12) 

and will not be repeated here. 

F. Improvement of Feasible Trajectory 
Equations (3.14) to (3,17) used to reconverge the trajectory at 
each iteration of the improvement phase must also be modified because 
of the variable model switching times. Suppose there is a model switch- 
ing time such that t^^ < t^. Then, if t^^ is modified, the effect 
on the model switching function S is given by 


AS = a-^ Ax(t ) 
s 


(4.31) 


where 4x(tg) is obtained by integration of (3.16) from t^^ to t^, 
with initial conditions given by (3.14). The change in S must be 
canceled by a change in t^. The result is 


-a'^ Ax(t^) 


(4.32) 


a x(t ) 


The terminal conditions are altered due to the direct effect of 

the change in t^^, and indirectly due to the change in t^. The 


total change in is given by 




Aijj, = 


9x^ a Ax(tg) 

Ax(tp - 3^ T...-V 

s a x(t ) 

s _J 


(4.33) 


where Ax(tp is obtained by integration of (3.16) from t^^ to t^, 
with initial conditions given by (3.14) , and 9x^/9tg is obtained by 
integrating (3.1a) subject to initial conditions (3.1b), for t^ = t^. 

The change in must be canceled by appropriate changes in 

the t^. As in (3.17), this results in 



CHAPTER V. DEMONSTRATIVE EXAMPLE 


The ideas developed in chapters II through IV are illustrated in 
this chapter with the aid of numerical examples. First, a problem 
will be solved in which the system equations are nonlinear in the 
state and linear in the control. The possibility of minimizing singu- 
lar arcs will be discussed. Next, the same problem will be solved by 
using one- and two-piece linear approximations to the original system 
equations. Solutions will be obtained for two different sets of 
terminal state constraints , and both one and two control variable 
problems will be solved. 

A. Problem Statement 

Consider the problem of finding u(t) which transfers the system 


X 



+ u 


( 5 . 1 ) 


y = -4y + u 

subject to the control limits juj —2 from initial conditions 
(xQ,yQ> = (1,1) to terminal conditions (x^,yp = (0,0) in minimum 
time. 


We have in (2.1), (2.3), and (2.4) 

(j) = t^, a = b = 0 

G- = -2 , d^ = -1 

(5.2) 

^^2 d ^ = 1 

” ^f ’ ^^2 ~ ^f 
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The Hamiltonian (2„5) is 


H = X (-x^ + u) -f > (-4y + u) + y, (-2 - 


and the Euler Lagrange equations are 


X = 2X. X 

X X 


X. =• 4X. 

y y 


u) + + u) (5,3) 


The control is determined from 


min(X + X )u 
X y 
u 


(5.4) 


(5.5) 


which results in 


u =■ -2 sgn(X + X ) 
X y 


The terminal conditions on the Lagrange multipliers are 

X(tp = V 


HCt.) = 


(5.6) 


(5.7) 


B. Nonexistence of Singular Arc 

The optimal control Is determined from (5.6) unless (X^ + X^) s 0. 
In this case, the control is determined by successively differentiating 
(X + X ) with respect to time until the control appears explicitly. 

X y 

This procedure results in 


4 - + A > - 2X X + 4X - 0 

dc - X V X y 


(5.8) 


(X + X - 2X x^ + 2X u -h 16X =0 

^j.2 X y X X y 

Simultaneous solution of ( + ■'^) — 0 eed (5,8) results in two possi- 
bilities for a singular arc. First, we may have x = 2 and u - 4. 
However, this is not possible, since we must have |uj _ 2, The second 


possibility is that X = X = 0. In this case, we would have H = 0 

X y 

and since the system is autonomous, H = 0. However, this is incon- 
sistent with H(tj) = -1 in (5.7). Therefore, there can be no singu- 
lar arc for this problem. 


C. Initial Feasible Solution 

Since there are two terminal constraints and the terminal time 
is free, the initial feasible trajectory must have two segments. The 
only possible control for these segments is u = ±2. Therefore, we 
assume tentatively that the optimal control history is given by 

u = -2 , 0 £ t £ t^ 


(5.9) 


u = 2, t^ £ t £ t^ 

with and t^ to be determined such that = y^ = 0. The state 

equations (5.1) can be integrated in closed form when (5.9) is used 


as the control. The result is 


1-/2 tan /2 t- 

x(t^) = ^ _■ 

1 + ™ tan t, 

/z ^ 


y(t,) = - I (l - e 


Iterative solution of (5.10) for t^ and t^ 


4 - 

+ e 

such that 


x(tp = y(tp = 0 results in 


0.56, 


t^ = 0.69 


(5.10) 


D. Calculation of Lagrange Multipliers 
•For this problem, we have 9(j>/9x = 0 and 9i|^/9x — I. Therefore, 
X(t ) = V. We must integrate (5.4) backward with three sets of initial 


conditions - X 


( 0 ) 


(tp 


0 
0 


(tf) 




in order to 


Equations (5 a 4) may be integrated in 


find A(t^) and 
closed form as follows : 




A 




-2X (t ) - 2A^(tJ - X (t)[-4y(t) + 2] 
X„(t) = ^ 


-x"(t) + 2 


t Z. t. 


Therefore, we obtain 




0 /l.O35\ 




" ) 

,0.595/ 


Equations (3.8) and (3.9) give 


4 (1,1) 


1.035 


0 N (^A 

0.595/ \v„/ 


= 0 


and 


( 2 , 2 ) 


= -1 


which can be expanded to give 


1„035 + 0.595 ^2 " ° 


Vi + V£ = -0,5 


Solving for and ^2 yields 

= 0,676, V 2 = -1.176 

Therefore, from (5.7) we have 

X (t.) = 0.676, X Yt ) = -1.176 (5.11) 

X f y t 

From this problem, the state equations and Euler-Lagrange equations 
can be integrated in closed form. The result is 


1 


F 


49 


Ay(t) = -1.176 e 


-1 - ^^(t)[-4y(t) + 2] 

A (t) = 2 

X 2, s 

-X (t) + 2 

-0.252 + 1.414 tanh (t - t ) 
x(t) = 1 

1 “ 0,178 tanh /2 (t - t^) 


\ < t < (5.12) 


r -Ht-tJ] 

y(t) = 0,5 1 - e - 


0.340 e 


-4(t-t^) 


Ay(t) = -0.699 e 


ACt-tJ 




1 - ^^(t)[4y(t) + 2] 
A,(t) V- 


X 


x^(t) + 2 




x(t) = X - 1.414 tan /2 t 
1 + 0.707 tan /2 t 

y(t) = -0.5(1 - e“^^) + 


(5.13) 


E. Optimality of Initial Feasible Solution 
The initial feasible trajectory was obtained under the tentative 
assumption that the optimal u(t) is given by (5.9). According to 
(5.6) y this control strategy is optimal only if 


A + X > 0. 
X y ’ 


A + A < 0 , 
X y 


0 <. t 


t ^ t t 

L-i _ c _ 


(5.14) 


To determine if this is the case^ (A + A ) is calculated as a function 

X y 

of time from (5.12) and (5.13), The result is plotted in figure 5.1, 
and shows that for this problem, the necessary conditions for opti- 




nL i 
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mality are satisfied by the initial feasible trajectory. Therefore, 
no further improvement needs to be made. The optimal trajectory, 

x(t) and y(t), is shown in figure 5.2. 

F. Linear Approximation 

The nonlinear system equations (5.1) may be approximated by 

linear equations. Actually, since the y equation is already linear, 

only the x equation needs to be approximated. If we choose to 

linearize the x equation about x = 0,5, the result is 

X = -2x| (x - 0.5) + (u - 0.25) 

|x=0.5 

= -X + u + 0.25 C5.15) 


For this linearized model, the Hamiltonian is 

H X (-X + u + 0.25) + X (-4y + u) + v- (-2 - u) + (5d6) 

X y ^ 

and the Euler-Lagrange equations are 



(5.17) 


As previously, the initial feasible trajectory is obtained by assuming 
(5.9) is the optimal control strategy. Solution of the resulting 
boundary value problem for t^^ and t^ yields 

t^ =0.68, t^ = 0.827 

The solution of equation (3.7) is given by 

/ 0.863 0 \ 

A(t.) = ( 

^ \0 


0.555 




and equations (3.8) and (3.9) are evaluated to give 

= 1.163, ■''2 ” -1.808 

from which it follows that 

X^(tp = 1.163, Xy(tp = -1.808 

Since the boundary conditions are known, the system and Euler-Lagrange 
equations can be integrated to give x, y, X^, and X^ as functions 
of time. The switching function (X^ + X^) is shown in figure 5.1, and 
the optimal trajectory is shown in figure 5,2. As was the case when 
the nonlinear equations were usedj, the initial feasible trajectory is 
a local minimum in the present case also. 

G. Piecewise^-Linear Approximation 
Comparison of the nonlinear and approximate linear results in 
figures 5,1 and 5.2 shows that the linear approximation does not give 
a very accurate representation of the nonlinear system for the present 
problemo The approximation may be improved by using a piecewise 
linear model,, as described in chapter IV. If we choose x - 0,25 and 
X = 0,75 as the equilibrium values about which linear models will be 
derived, the resulting linear models are 

X = - y + u + Yg- ^valid near x = 

(5.18) 

X, = - Y X + u + ^ ^valid near 

Since the y equation is linear, the weighting matrix W in (4,15) 


is chosen to be 
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This choice of W results in a model switching function (see (4a6) 
and (4.17)) of 


S = X 


(5 c 20) 


The Euler-Lagrange equations are 

' I 


X 1 0.5 


^x - 2 ^x> 


X ^ 0.5 


(5.21) 


.-jr'*-. 


X = 4X 

y y 

At the model switching point, there is a jump discontinuity in 




( 5 . 22 ) 


where e is determined from (4,20) 


e = 


X.( Of-7 (0-» -1J. . 

i(0.5) +2+^ 


Which, results in e - 0 for this problem. 

Based on the earlier solutions presented in figures 5.1 and 5.2, 
it appears that the model switch will take place prior to the change 
of control at t = t^. Therefore, we assume that t^ < t^ and look 
for a feasible solution with control given by (5.9). Solution of the 

multipoint boundary value problem gives 

t = 0.197, t-. = 0.597, t^ = 0.733 

s 

The solution of equation (3.7) is 


0.9H 0 


A(tf) 
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and equations (3o8) and (3o9) are solved to give 


0,863, ly(tp = V2 = -1.390 

Again, the trajectory and Euler-Lagrange equations are solved as 
a function of time, and the results are presented in figures (5.1) 
and (5.2), It can be seen that the two-piece linear model is a good 
representation of the actual nonlinear system for the present problem. 

H. Additional Control Variable 
Consider the following problem. Find u(t) and v(t) which 
transfers the system given by 


x=-x+u+yv+0,25 

X 

(5.23) 

y = -4y + u + YyV 

from initial conditions (xq, Yq) = (1,1) to terminal conditions 
(Xf,yf) = (0,0) in minimum time. The controls are subject to the 
constraints 

juj i 2 and 0 V ^ 1 

This problem is identical to the one-piece linear approximation 

of the original nonlinear example problem, except that there are now 

two control variables, u and v. The coefficients Y and y will 

X y 

be left unspecified temporarily. 

Since there are two terminal constraints, the initial feasible 
trajectory must have two segments. The allowable control strategy for 
either of the two segments is u = +2 or -2 and v = 0 or 1< Sup- 
pose we choose for our control strategy 
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u = -2 and v = 0 , 0 £ t 1. 


(5.24) 


u = +2 and v = 0, 


+- ^ t" ^ t 


T, IT = 0 for both segments, the initial feasible 

jince we have chosen v 

n-n-T^ii<sr time histories are identical to 
trajectory and Lagrange multxplxer t 

T-iPr for the one-^iece linear model, single control 
those obtained earlxer for tne one . 

j. "x ^ for control n is shown in 

problem. The switching functxon 

figure 5.1, and the trajectory Is shown in figure 5.2. 

OTT, tjp are interested not only in the opti- 
For the present problem, we are 

mallty of uj^(t) , but also in the optimality of Vj^(t). The switc 

ing function for v(t) is given by + Y/,)1 optimal v(t) we 


nust have 


V = 0 if (Y^x^ + YyXy) > 0 

V - 1 If (Y,^, + YyXy) < 0 


(5.25) 


f \ and X are shown in figure 5.3. The xnx- 
Che time hxstorxes of and 

tial feasible trajectory is a local minimum if the values of Y^ 

T are such that Cy,X, f Y^Xy) > 0 for the entire trajectory; for 
instance, this is the case if 1 and - 0.5. 

hand, if Y^ = 1 and Yy = 0.8, the switching function (X^ + 0.8 Xy) 

is Shown as a function of time in figure 5.4. For this case, v = 0 

, . _ n 75 sec but the negative 

is the optimal control up to about t - . 

< £ ... f- t> n 75 sec indicates that t^ 

value of the switching functxon for t * 

can be reduced if v is increased in this re gxon. 

in order to improve the trajectory (decrease the performance 
index, t,) a third control se^ent is added to the trajectory. The 
rh;rd secment is «de very short initially, so that the new trajec- 
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toTy does not differ eubstantielly from the initial feasible trajeo- 
tory. Specifically, the control is chosen to be 


u = “2 and 

V = 0, 

0 .1 t i t 

u = +2 and 

< 

II 

o 

ti < t £ t^ 

u = +2 and 

< 

II 

1-* 

t 1 t 1 1 



sw 

is selected. 

and held 

fixed while 


varied to reoonverge on the end conditions . 0. The converged 

values of and t^ are given by 

tx * On 68, t^ = 0.819 

It can be seen that t^ has been reduced from its value of 0.827 for 
the initial feasible trajectory. At this point, the Lagrange multi- 
plier time history must be calculated for the new three segment tra- 
jectory, and the gradient-of the performance index with respect to 
'^sw dalculated by using (3.18). A nonlinear search technique such as 
presented in reference 11 is used to iteratively search for the value 
sw which t^ is a minimum. When this has been accomplished, 

the switching functions (1^ + and (1^ + 0.8 X^) should be exam- 
ined once again to determine if t^ can be further reduced by includ- 
ing additxonal control segments in the trajectory. 



I 


1 


¥ 

4 


I 


I 


CHAPTER Vic APPLICA.TION TO TURBOFAN CONTROL 
In this chapter, the algorithm described in chapters III and IV 
is used to find optimal trajectories for an FlOO aircraft turbo fan 
engine o Specifically, values of the control variables are found as a 
function of time, which minimize the terminal error for a fixed-time 
acceleration, while adhering to the engine constraints. This is equiv- 
alent to minimizing acceleration time for a fixed terminal error, A 
suboptimal, closed loop control mode is also developed. Finally, the 
problem of minimizing steady-state specific fuel consumption is con- 
sidered, 

A. Engine Description 

The following description of the FlOO engine is taken largely 
from reference 18, The Pratt & Whitney FlOO engine (fig. 6.1) is an 
axial, mixed-flow, augmented, twin-spool, low-bypass-ratio turbofan. 

A single inlet is used for both the fan airflow and engine core air- 
flow. Airflow leaving the fan is separated into two streams : one 
stream passing through the engine core and the other stream passing 
through the annular fan duct, A three-stage fan is connected by a 
through-shaft to the two-stage, low-pressure turbine. A ten-stage 
compressor is connected by a hollow shaft to the two-stage, high- 
pressure turbine. The fan has variable, trailing-edge inlet guide 
vanes, whxch are positioned by the engine control system as a function 
of fan corrected speed to maintain fan stability at low speeds. The 
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rHIGH-PRESSURE TURBINE 



figure 6.1. - Schematic representation of ROO-PW-lOO augmented turbofan engine. 
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compressor has a variable inlet guide vane followed by two variable 
stator vanes; the vanes are positioned as a function of compressor cor- 
rected speed. The engine core and fan duct streams combine in an aug~ 
mentor and are discharged through a variable convergent-divergent 
nozzle. 

The fuel control system meters fuel to the main combustor, as a 
function of power lever angle (PLA) , compressor speed, fan discharge 
temperature and compressor discharge static pressure. Augmentor fuel 
flow is metered as a function of PLA, fan discharge temperature, and 
compressor discharge static pressure. Exhaust nozzle area is con- 
trolled so as to maintain a desired engine airflow during augmented 
opei:v^.tiono 

B. Engine Models 

Pratt & Whitney Aircraft (P&WA) has developed a detailed dynamic 
simulation of the FlOO engine using a digital computer (ref. 14). 

The simulation xncludes overall performance maps of the engine compo- 
nents, Variable gas properties, and Reynolds number effects in order 
to provide good steady-state accuracy over the range of power settings 
and flight conditions. Factors such as fluid momentum, mass and energy 
storage and rotor inertias are included to provide transient capability. 
A detailed simulation of the engine's control system is also included. 

In addition to transient capability, the simulation also has the capa- 
bility to solve iteratively for ecjuilibrium operating points for spec- 
ified flight conditions and power lever angle, 

A drawback to the use of detailed dynamic simulations on a digi- 
tal computer is that they require large amounts of computer time to 
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obtain transient solutions. A requirement of several minutes of digi- 
tal computer time per second of real time is typicals Szuch and 
Seldner (ref. 18) developed a hybrid computer simulation of the FlOO 
engine which runs in real time. This simulation is useful for develop- 
ment and checkout of engine control system hardware and software. 

Some simplification is made in the hybrid model 'in order to allow the 
real time capability. Nevertheless, steady-state and transient re- 
sults presented in reference 18 compare favorably both with the de- 
tailed P&WA digital simulation and with a limited amount of experi- 
mental data. 

Because of the virtual impossibility of using nonlinear feedback 
control theory for realistic systems, control software is usually 
developed using linear models. For turbofan control system design, 
nonlinear dynamic simulations such as references 14 and 18 are linear- 
ized about various equilibrium conditions, and linear models obtained. 

This process produces equations of the form 

X - F(x x^) + G(u - u^) (6*1) 

where x^ and u^ are equilibrium values of state and control, re- 
spectively. Other engine variables which are not modeled as states 
are also of interest. Such variables will be called outputs, and 
denoted by y. The linearized output equations are given by 

y = y + C(x - X ) + D(u - u ) (6«2) 

0 C 6. 

If the Outputs have upper (or lower) bounds which must not be exceeded, 
then combined state/control path inequality constraints of the form 
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+ C(x - x^) + D(u - u^) - i 0 

(6.3) 

Vn • 5'e - ' “e> - “ 

are produced. Mechanical limits on the control variables also have 
the general form of (6.3) with C = 0. 

The engine simulations of references 14 and 18 contain 16 and 
17 state variables, respectively. Therefore, the state vector x in 
(6.1) has dimension 16 or 17, depending on which simulation is used 
to obtain the linear model. It is preferable to conduct control sys- 
tem studies by using lower-order models, if possible. Fortunately, 
in the present case, most of the eigenvalues (natural frequencies) of 
the system matrix F are considerably larger than the lowest eigen- 
values. Therefore, quasi-steady approximation may be used to reduce 
the order of the system. A quasi-steady approximation technique which 
preserves the lower eigenvalues exactly, and also retains the desired 
states, is presented in appendix C. 

In a recent contractual effort under joint Air Force/NASA sponsor- 
ship (ref. 8), Systems Control Inc. (SCI) used linear quadratic regu- 
lator theory to design controls for the FIDO engine. Linear models 
having 16 states and reduced models having 5 states were provided to 
SCI by P&WA for a number of equilibrium points at different flight 
conditions and PLA’s. Some of these linear models are given in refer- 
ence 13. 

Five of the five-state equilibrium models from the P&WA/ SCI study, 
equally spaced along the sea-level-static (SLS) operating line, are 
used in the present report. The normalized linear model for 
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PLA = 67 degrees is shown in table I. The 67 degree power lever set- 
ting is typical of subsonic cruise; PLA is 20 degrees at idle, and 
83 degrees at maximum nonaugmented thrust (also referred to as inter- 
mediate) ■> The normalized equilibrium values of the states, controls, 
and outputs are given in table II, along with the definitions of these 
variables . 
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TABLE I. - FIVE-STATE LINEAR MODEL 
[Sea level static, PLA = 67 °] 


-4,20 

0.877 

-0 .645 

1.22 

-0.353 

-5.63 

-0.335 

-0.416 

19,2 

-13.1 

-7.38 

-0.257 

6.28 

-21.9 

-3.02 

-30.2 

71.0 

456 

39.2 

-1.06 

-0.573 

-0.0349 

-0 . 103 

-0.000639 

.976 

-0.0212 

-0,0256 

-0.00591 

2,11 

-7.50 

.702 

.0138 

32 

-0.281 

-0.174 

.0351 


- 0.666 

0.0292 

-0.349 

2,31 

23.7 

-0,0135 

.0432 

-0.468E-2 

.139 

.327 


0.154 

-0.0194 

-0,0498 

-1.38 

.776 

0,0294 
.0799 
-0.30 3 E-2 
-0,669 
-0.0187 


0.109 / 

-0.0335 

0.196E-3 

.527E-3 

.0148 

.874E-2 

-0.614E-4' 

-0.560E-3 

.392E-3 

-0.438E-2 

-0.107 


0.892 

.813 

4.74 

3.22 

■147 


0.222 > 
-0.00131 
.0521 
-0.0963 

-3.62 i 


In order to determine if the fifth-order models could be fur- 
ther reduced, the eigenvalues and eigenvectors of the fifth— order 
models were calculat&d.. The results are presented in table III for 
PLA =67 degrees. The eigenvectors are the columns of the symmetry 
transformation T as defined in appendix C. It can be seen that there 
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TABLE II. - EQUILIBRIUM VALUES OF PROBLEM VARIABLES 


[Saa level static, PLA = 67°] 

Fan speed, 100 rpm \ 

Compressor speed, N , 100 rpm 

comp ’ ^ 

Augmentor pressure, Pt7, 0.1 psi 

Fan turbine inlet temperature, FTIT, 10'’ R 

Combustor pressure, P^^ » J 

Combustor fuel flox^, W^, 100 lb m/hr 

Exhaust nozzle area, A , 0.01 ft^ 

noz 

Inlet guide vane position, IGV, 0.1 deg 


High variable stator position, HVS, 0,01 deg 


94,39 


121. 70 


187.1 


68.60 


Thrust, T, 100 lb 

Airflow, w^, Ibm/sec 

Turbine inlet temperature, TIT, 100% 

Fan surge margin, SMFAN (0.001) 

Compressor surge margin, SMCOMP (0.001) 


105.14 


86.42 


TABLE III. - EIGENVALUES AND EIGENVECTORS OF FIVE-STATE MODEL 
[Eigenvalues = -3.42, -4.42 # 2.91j, -30.8, -151 (sec“^) ] 


T 


-0.0823 0.317±0j -0.0463 0.00594 

-0.118 0.0508+0. 0117 j .0152 .00574 

-0.761 1,0+1. 8j .0496 .0326 

.0884 0.00625+0. 142 j .997 .0281 

-0.626 -0.587+0. 470j .0392 -0.999 


-0.999 
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are three dominant eigenvalues, one real and a complex pair; the other 
two eigenvalues are at least a factor of six larger. Also, it appears 
from the T matrix that states four and five (FTIT and Ps3) are 
closely related to the two larger eigenvalues. Therefore, quasi-steady 
reduction was performed on the five state system, as described in 
appendix C, to reduce it to a three state system, with state variables 
^fan* ^comp* Pt7. It should be noted that these are the same 

three states used by Weinberg (ref. 6). The reduced three-state linear 
system matrices are presented in table IV. There are seven outputs for 


TABLE IV. - THREE-STATE LINEAR MODEL 
[Sea-level static, PLA - 67 °] 


F = 


3.13 

-2o87 

1.93 


-0.493' 

-0.072 

-6.03 


X 



~ 1 ^''comp 



G = 


C = 


D = 


-0.0354 

.0165 

-7.34 

0.0511 

.0249 

-0.187 

2.00 

12.3 

-0.465 

3.16 

-0.00523 

.0432 

-0.00294 

.135 

.206 

-0.025 

.0335 


-0.0848 

-0.00795 

.804 



^0.00155 

-0.00574 

.0173 


ye 


\ 


^a 
TIT 

SMFAN 
SMCOMP 
FTIT 
,PS3 


1 


0.0342 

.0799 

-0.00191 

-0.671 

-0.0967 

™0.437E~3 

.0215 


0.807E-4 
-0.561E-3 
.436E-3 
-0.444E-2 
- 0.110 
0.119E-2 
.819E 


-I! 
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the three state model, since the variables FTIT and Ps3 now have the 
role of outputs. The eigensystems for the four other linear models 
have the same qualitative features as observed in table IV, and the 
same quasi-steady reduction was perforiEed on these models also. 

C. Model Accuracy 

A nonlinear model can be approximated by a suitable number of 
linear models at various equilibrium points. It is of interest to 
determine how many equilibrium models need to be used for suitable 
accuracy. One way to do this is to see how well the equilibrium con- 
ditions for one model are. predicted by adjacent models. For example, 
if the equilibrium control from model 2 is used in model 1, the pre- 

/s 

dieted equilibrium state for model 2, ±s given by 

which implies 


(6.4) 


X r, = X 1 

e2 el 


The error in predicted state is 


F^ G^(u,o - 


-1. 


(6.5) 


= ^e2 ■ ^e2 == ^^el ' ^e2^ h 


( 6 . 6 ) 


The error in predicted output can be calculated by using ^^2 

X ,, from (6.5) in (6.3), The result is 
e2 ' 




'e2 ■ 'el ' "“1 

The error in predicted output is 

~ ^e2 " ^ye2 " ^el^ ^^1 " Vl^^l^ ^%2 “ 

The results obtained by using (6.4) and (6.8) with the five 
linear models are presented in figure 6.2. The true values of the 
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problem variables (states and outputs) are presented as a function of 
PLAo Tlae values at the five PLA’s are connected by straight lines; 
however, the true variation between equilibrium data points is not 
necessarily linearo The predicted problem variables are also shown, 
as predicted from adjacent equilibrium points both above and below 
the predicted pointy 

It can be seen from figure 6,2 that prediction accuracy is gen- 
erally very good, with, several notable exceptions. First, for nearly 
all of the variables, the idle model (PLA = 20 degrees) does not accu- 
rately predict the equilibrium values at PLA = 36 degrees. Also, the 
model at PLA = 36 degrees does not predict the idle conditions accu- 
ratelyo Both of these results may be explained by the fact that the 
engine exhaust nozzle and low pressure turbine are unchoked at idle, 
but become choked a few degrees above idle and remain choked as PLA 
increases further. Engine dynamic characteristics differ substantially 
between choked and unchoked conditions. Therefore, there is a large 
discontinuity in the linear models at the point at which choking occurs. 

With the exception of the idle prediction difficulty, all problem 
variables are accurately predicted over the range of PLA, except for 
the fan and compressor surge mar gins. The difficulty in predicting 
surge margins is due to the fact that surge margin is proportional to 
the difference between two pressures which are of similar magnitudes. 
Hence, relatively small errors in the pressures give rise to large 
errors in the surge margins. 

Because of the fact that the engine exhaust nozzle is unchoked 
only in the near vicinity of idle, and furthermore is choked at all 
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PLA's at nonzero flight Mach numbers, it was decided to eliminate the 
idle model from all subsequent calculations. Therefore, the engine 
will be modeled by using the three— state models (e.g, , Table IV) at 
four equilibrium points: PLA == 36, 52, 67, and 83 degrees. 

D. Steady State Performance 

Linear models can be used to find the control settings which 
maximize steady-state performance. For example, suppose we wish to 
maximize thrust for constant fuel flow — this is equivalent to minimiz- 
ing specific fuel consumption. The minimization must be accomplished 
while adhering to the engine constraints and control limits. The mini- 


mization will be conducted at PLA - 67 degrees, a typical power setting 
for subsonic cruise. At this power setting, the only applicable engine 
constraints are that fan and compressor surge margins should be kept 
above a safe level, say five percent. In addition, the exhaust nozzle 
area, inlet guide vanes, and compressor variable vanes must not exceed 
their mechanical limits , 

For given values of u, the equilibrium values of states and out- 
puts are given by (6.5) and (6.7) 


X = X - F ^G(u - u ) 
ss e e 


y ” y 

''ss e 


+ (D - CF“"'G) (u - u^) 


(6.9) 


If the performance "'ndex and constraints are linear combinations of 
States, outputs, and controls, a linear programming problem results, 
which can be solved by usxng the simplex method (ref, 16) . 

Results for the maximization of thrust are presented in table V. 
For all cases, fuel flow is held constant at 6860 Ibm/hr, For the 
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TABLE V. “ STEADY-STATE PERFORMANCE 


[Sea-level 

static. 

PLA = 

67% W£ 

= 6860 Ibm/hr] 



A , 

IGV, 

HVS 

Thrust , 

SMFAN 

SMCOMP 


noz ’ 
2 

deg 

deg 

Ibf 




ft 






Nominal 

2.98 

-16.4 

0.92 

10 514 

0.25 

0.18 

Optimum 

2„8 

5 

-40 

11 411 

.07 

.71 

Scheduled geometry 

00 

CM 

-9.4 

-6.1 

10 790 

.16 

00 


±7° 

nominal case, the exhaust nozzle area, inlet guide vanes and compressor 

variable vanes are at their scheduled values; the thrust which results 

is 10 514 Ibfo In the optimum case, A , IGV, and RVS are free for 

noz 

optimization, subject to their mechanical limits. The resulting 
thrust is 11 411 Ibf, which represents an increase of about 9 percent 
over the nominal value. However, it is recognized that the linear 
models are not necessarily valid over the full range of possible con- 
trols. Furthermore, use of optimum values of IGV or HVS could 
result in violation of flutter boundaries. Therefore, a third case is 
considered in v?hich the values of IGV and HVS are allowed to devi- 
ate by at most 7 degrees from the scheduled values at that power set- 
ting. In this case, the thrust is 10 790 Ibf, which still represents 
an increase of more than 2.5 percent from the nominal value. 

A possible explanation for the increased thrust achieved by using 
off-schedule geometry is that the scheduled values of IGV and HVS 
are based on optimized component performance, rather than optimized 
overall engine performance. 
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E. Transient Performance 

We now consider the problem of minimizing the time required to 
accelerate the FlOO engine from one equilibrium thrust level to another 
as rapidly as possible. In solving this problem, a piecewise-llnear 
model of the FlOO engine will be used. Specifically, the engine model 

is given by 

X = F^(x - + G. (u - (6.10) 

where the state vector x (three states) and control vector u (four 
controls) are as defined in table IV. There are also other problem 
variables of interest, called outputs, which are related to the states 

and controls by 

y - C^(x - X^.) + dJ(u - U^,) + y^^ (6.11) 

the index 1 refers to the equilibrium model number; there are tour 
equilibrium linear models at FLA - 36, 52, 67, and 83 degrees. The 
67 degree model is presented in table IV. The linear model which 
applies at a given time is selected by minimizing the quadratic func- 

tion 

= (X - x^i>^W(x - x^i) (6.12) 

with respect to i; the model whose equilibrium state is "closest” to 
the actual state at that time is chosen to represent the engine. The 
function of W in (6.12) is to scale the states to comparable numeri- 
cal values 0 

During an acceleration from a part-power condition to intermediate, 
the engine is first represented by the PLA = 36 degrees model. As the 
engine accelerates , the model switches successively to the 52, 67, and 


83 degree models. The W matrix used in comparing successive PLA 
models k and (k + 1) in (6.12) is given by 



It should be noted that the components of W defined by (6.13) are 
different for each successive switch. The normalizing factor for 
each state has been chosen to be the difference in values of that 
state from one equilibrium point to the next. 

The trajectory must also satisfy path inequality constraints given 
by 

cTx + d^u + e,<.0, i = l, . . .,q (6.14) 

XXI 

The coefficients c,, d. , and e. are different for each equilibrium 

X , x’ X 

model. Some of the path inequality constraints correspond to engine 
physical limits, others to control mechanical limits. The following 
constraints will be assumed for this problem. 

(1) Turbine inlet temperature cannot exceed the equilibrium value 
at intermediate thrust by more than 50 degrees !R. 

(2) Fan and compressor speeds cannot exceed the equilibrium values 
at intermediate thrust by more than 50 rpm. 

(3) Fan and compressor surge margins must not be less than 5 per- 
cent. 

(4) Inlet guide vanes, compressor vanes, exhaust nozzle area, and 
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fuel flow rate must not exceed their limits. 

Statement of the problem . - Before developing a precise mathe- 
matical statement of the optimization problem to be solved, it is use- 
ful to consider first the manner in which the resulting optimal control 
for minimum-time acceleration (which will be referred to as the transi- 
tion control) might be implemented in an actual engine. The transition 
control will be used during acceleration from one equilibrium power 
setting to another; during near-steady conditions, the Bill-of-Material 

ic 

(BOM) control (which will be referred to as the regulator) will be 
used. 

In deciding when to transfer authority from the transition con- 
troller to the regulator, it must be recognized that the equilibrium 
conditions at a given PLA are not known precisely ; they vary from 
one engine to another and as a function of operating time for a given 
engine. Therefore, control authority should be transferred from the 
transition controller to the regulator when the state vector is in the 
"vicinity” of the desired state vector rather than when desired values 
of the states have been achieved precisely. The distance to the de- 
sired state vector (terminal error) may be defined as 


( 6 . 15 ) 


The switch from the transition conttoller to the regulator would occur 



The BOM control is the standard control system supplied with the 
engine. 
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vj ) at a fixed value of say 0.05 or 0.1 

The problem to be solved is stated as follows. Find controls 
u(t>, 0 <. t 1 t^ which minimize <f,(xp while satisfying the system 
equations (6„10) and path constraints (6«14). A sequence of solutions 
to such problems for different acceleration times may be used to find 
the minimum-time solution for a given value of terminal error. Neces- 
sary conditions for an optimal solution are given in chapter II. The 
problem is solved by using the new algorithm described in chapters III 
and IVc 

Results o - The problem of minimizing the terminal error for an 

acceleration from near-idle (PLA =24 degrees) to intermediate thrust 

is considered. The engine's exhaust nozzle first becomes choked at 

PLA = 24 degrees. The final time is specified to be t^ = 0.75 second. 

The problem variables (states, outputs, and controls) for the optimal 

trajectory are shown in figure 6.3. The state variables N N 

fan * comp ’ 

and Pt7, are shown as functions of time in figure 6.3(a) to (c) . It 

can be seen that the states approach the desired final values smoothly 

and with no overshoot. The value of the terminal error is ^ = 0.012, 

which results from errors in the states of AN = 15 mm 

fan ^ » 

The outputs are shown as functions of time in figures 6.3(d) to 
(j)o Because of the way in which the outputs are defined in equa- 
tion (6„11) these variables are in general discontinuous at model 
switching points and points of discontinuous control. 

The optimal control strategy results in the high-pressure turbine 
temperature having its maximum value for the entire trajectory^ 
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this is shown in figure 6.3(d). Fan and compressor surge margins 
remain well above acceptable minimums (figs, 6.3(e) and (f))t Low 
pressure turbine inlet temperature (fig. 6.3(g)) is very nearly con- 
3 tant. Thrust (fig. 6.3(h)) increases smoothly and laonotonlcally 
with time. 

The optimal control histories are shown as a function of time in 
figures 6.3(k) to (n) . Fuel flow jumps at t = 0 from its idle 
value, then increases slowly to maintain constant turbine inlet tem- 
perature. The optimal values of piecewise 

constant, as required by variational theory (chapter II). Each of 
these variables has one switch during the trajectory. 

Figure 6.3(o) shows the distance from the current state vector 
to the equilibrium state vector of the current model. The distance 
is normalized in such a way that the distance between adjacent equi- 
librium state vectors is unity. It can be seen that the instantaneous 
distance is always less than unity; this is a good indication of the 
validity of the piecewise-linear model throughout the entire trajectory. 

Although the maximum nozzle area is greater than 6 square feet, 
nozzle area has been restricted to a range of 2.8 to 3.2 square feet 
for the present study. This arbitrary upper limit has been imposed 
because the linear models are not necessarily valid for the full 
range of allowable control action. It is also recognxzed that model 
accuracy may be degraded if IGV and HVS values are far. from their 
scheduled values from the BOM controller. In fact, available test 
data for the. IGV and HVS are limited to about ±7 degrees of the sched- 
uled values. Furthermore, large deviations in IGV and HVS might result 
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in flutter "boundaries being violated. For these reasons, two addi- 
tional optimal trajectories were run for t^ = 0.75 second. For one 
of the trajectories, the HVS and IGV were required to have their 
scheduled values; for the other, IGV and HVS were limited to ±7 degrees 
from their scheduled values. In addition, optimal trajectories were 
run for other values of t^. 

Figure 6.4 shows the problem variables as a function of time for 
the optimal trajectory in which the IGV and HVS were limited to ±7 de- 
grees of the scheduled values. The acceleration time is 0.75 second. 

For this case, the terminal error is 0.045, compared to a terminal 
error of 0,012 which is achieved when the IGV and HVS are allowed full 

yaxiation wi.thin their mechanical limits . 

It can be seen that the results for this case are qualitatively 
very similar to those of figure 6.3. Turbine inlet temperature has 
its maximum value for the entire trajectory, and the state variables 
increase monotonically to their final values, IGV position 
(fig. 6,4(m)) is 7 degrees less than the scheduled value up to about 
0.56 second; after that time, IGV position is 7 degrees greater than 
the scheduled value. It is interesting to note that HVS position 
(fig. 6.4(n)) is 7 degrees less than the scheduled value for the 

entire trajectory. 

Figure 6.5 presents terminal error as a function of acceleration 
time, t^, for accelerations from PLA = 24 degrees to intermediate. Re 
suits are presented for values of IGV and HVS which are fully opti- 
mized, scheduled and scheduled ±7 degrees. The corresponding trajec- 
tories have the same characteristics ps shown in figures 6.3 and 6.4. 
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(n) COMPRESSOR VARIABLE VANE POSITION 
Figure 6.4, - Concluded. 



TERMINAL ERROR 



Figure 6.5. ^ Terminal error versus 
acceleration time. Part-power 
(PLA = 24 deg) to intermediate 
thrust acceleration. 
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In particular, turbine inlet temperature is kept at its maximum value 
throughout each of the trajectories. 

As discussed earlier, authority might be transferred from the 
transition controller to the regulator at a fixed value of terminal 
error. The time required to reach a given terminal error depends on 
the control strategy used. Tor example, figure 6.5 shows that a time 
of 0.80 second is required to reduce the terminal error to 0,05 if 
scheduled IGV and HVS are used. If IGV and HVS are controlled opti- 
mally within scheduled ±7 degree values, the time is reduced to 0,74 
second; if fully optimized IGV and HVS are used, t^ is 0.65 second. 

The time required to accelerate from other initial values of PLA 
to intermediate is also of interest. Figure 6.6 presents the time 
required to accelerate for a terminal error of 0.05 as a function of 
PLA. Results are presented for scheduled and scheduled ±7 degree 
values of IGV and HVS. 

F. Comparison of Nonlinear and Piecewise-Linear Responses 
In figure 6,2, accuracy of the piecewise -linear model was investi- 
gated in terms of the ability of the model to predict the values of 
adjacent equilibrium values of states and outputs. This type of accu- 
racy test is particularly relevant to the ability of the piecewise- 
linear model to predict values of the control variables for optinial 
steady— state performance. However, for transient performance predic 
tion, a much better check on model accuracy is achieved by comparison 
of the nonlinear and pie cewis e-line ar model transient responses for 
the same control variable time histories. 

Figure 6.7 presents a comparison of nonlinear and piecewise-linear 
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Figure 6.6. - Time required to accelerate to intermediate thrust 
with a fixed terminal error of 0. 05. 















(i) FAN TURBINE INLET TEMPERATURE. 



(k) COMBUSTOR FUEL FLOW. 



(I) nozzle AREA. 
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transient responses. The control variable histories are based on an 
optimal acceleration from PLA = 24 degrees to intermediate thrust, 
for a specified acceleration time of 0.8 second. For this case, the 
values of IGV and HVS were constrained to be within +7 degrees of the 
BOM scheduled values. The nonlinear responses were obtained by using 
the P&WA digital dynamic nonlinear FlOO engine simulation (ref. 14). 

It can be seen that the nonlinear and piecewise-linear responses of 
compressor speed and augmentor pressure are in good agreement. How- 
ever, differences are observed in the fan speed responses. Also, there 
are substantial and important differences in the high-pressure turbine 
inlet temperature and fan and compressor surge margin responses. The 
nonlinear results show that the maximum value of the high-pressure 
turbine inlet temperature is violated by a large amount , and the com- 
pressor surges at about 0.06 second. The fan does not surge, but the 
fan surge margin does fall below the minimum value of 5 percent early 
in the trajectory. 

Additional comparisons of nonlinear and piecewise-linear results 
are made in figures 6.8 and 6.9. In these two figures, the comparison 
is based on an acceleration of 0.6 second from PLA = 36 degrees to 
intermediate thrust. For figure 6,8, the IGV and HVS are limited to 
scheduled values ±7 degrees, while for figure 6,9, the IGV and HVS 
take on scheduled values . The results are qualitatively similar to 
those presented in figure 6.7, but the nonlinear and piecewise linear 
results do not differ by as much as was observed in figure 6.7. For 
example, there is no compressor surge in figure 6.8, and neither fan 
nor compressor surge margin limits are violated in figure 6.9. 
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Figure 6.8. - Comparison of non- 
linear and piecewise linear re- 
sults. Acceleration from PLA = 
36 deg to intermediate thrust, 
acceleration time, 0. 6 sec. IGV 
and HVS, scheduled ±7 deg. 
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There are several possible explanations for the differences be- 
tween nonlinear and piecewise-linear responses observed in figures 6.7 
to 6.9. They are as follows: 

(1) The individual linear models may not be good approximations 
even for small perturbation inputs. 

(2) The number of equilibrium linear models may be inadequate to 
accurately describe the system nonlinearities with respect to the 
state variables . 

(3) There may be substantial nonlinearities with respect to the 
control variables , which are not included in the piecewise-linear 
model. 

(4) Model reduction from sixteenth to third order may have re- 
siilted in modeling inaccuracies. 

Items (1) and (4) may be checked by comparison of nonlinear and 
three and sixteen variable linear transient responses for small- 
perturbation control inputs. Such a comparison is made in figure 6,10, 
The equilibrium state and three and sixteen variable linear models cor- 
respond to PLA = 52 degrees, and the control input is a small step in 
combustor fuel flow. Results are presented for fan speed, compressor 
speed, augmentor pressure, turbine inlet temperature, and fan and com- 
pressor surge margins. It can be seen that the linear and nonlinear 
responses are in good agreement for states compressor speed and aug— 
mentor pressure, but are not in good agreement for fan speed. The 
sixteen-state linear result is in better agreement with the nonlinear 
result than is the three-state linear result, but neither of the linear 
results can be considered to be in adequate agreement. Also, the 
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linear and nonlinear results are in reasonably good agreement for tur- 
bine inlet temperature, but the three and sixteen state linear models 
do not give a good representation of fan or compressor surge margin. 

The differences between nonlinear and linear predictions of fan 
speed and fan and compressor surge margins observed in figure 6.10 are 
qualitatively similar to those observed in figures 6.7 through 6.9. 

It appears that a more accurate set of linear models, obtained from 
small perturbation responses of the nonlinear simulation, would result 
in more accurate piecewise -linear results. Also, there are substantial 
differences between the three and sixteen state linear results for sev- 
eral of the variables. The three-state linear models were obtained by 
modal reduction from five state linear models, which in turn were ob- 
tained directly from the nonlinear simulation. Since the three-state 
models were not obtained by modal reduction from the sixteen-state 
linear models, no conclusions can be drawn here relative to errors 
introduced by modal reduction. 

Because of the inaccuracy of the linear models, it is not possible 
to determine conclusively if nonlinearities with respect to states or 
controls contributed substantially to the differences in the nonlinear 
and piecewise -linear results observed in figures 6,7 to 6.9. Neverthe- 
less, the relatively good agreement of turbine inlet temperature ob- 
served in the small perturbation results (fig. 6.10) as compared to 
the large turbine inlet temperature errors observed in figures 6.7 to 
6.9 suggests strongly that there is a substantial nonlinearity, prob- 
ably with respect to combustor fuel flow, which would have to be in- 
cluded in the model to achieve accurate results. 


(b) COMPRESSOR SPEED. 
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TIME, sec 

(d) TURBINE INLET TEMPERATURE. 

Figure 6. 10. - Comparison of nonlinear and 
three-and-sixteen-state linear model 
responses to a three percent step in fuel 
flow. Equilibrium condition. PLA = 

CO J * 
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G. Comparison of Minimum Time and BOM Control Responses 
Using tlxe Nonlinear Simulation 

It is of interest to compare the transient responses obtained by 
using optimal minimum-time strategy with those obtained by using the 
BOM control. Such comparisons are made in figures 6.11 and 6.12, for 
accelerations from PLA = 24 degrees to intermediate thrust and from 
PL^ - 35 degrees to intermediate thrust, respectively. For both 
figures, the minimum time controls were obtained by using the piecewise- 
linear model with the IGV and HVS constrained to be within ±7 degrees 
of the BOM scheduled values. However, the data presented xn fxg- 
ures 6.11 and 6.12 were obtained by using the minimum-time and BOM 
controls applied to the nonlinear FlOO engine simulation. 

The results show that the minimum-time control strategy produces a 
more rapid acceleration to the vicinity of intermediate thrust. It 
appears that the principal reason for the improved acceleration is the 
much more rapid increase in fuel flow, which results in a rapid in- 
crease in high-pressure turbine inlet temperature. Naturally, the com- 
parison of performance is invalidated because of the violation of con- 
straints which occurs. Nevertheless, it seems highly probable that 
substantial improvement in acceleration time can be made without viola- 
tion of engine constraints since the high-pressure turbine inlet tem- 
perature (figs. 6 0 11 (d) and 6.12(d)) increases very slowly when the 
BOM control is used. Some compromise between the optimal control 
based on the piecewise-iinear model and the BOM control would probably 
yield improved accelerations without violating the turbine inlet tem- 
perature limit o 
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H. Suboptimal Control 

It has been shown that a substantial decrease in engine accelera- 
tion time can be achieved by using a minimum-time control strategy. 
However, the control strategy which has been developed is based on 
open-loop controlo It is desirable that a closed-loop control strat- 
egy be developed which is capable of closely approximating the open- 
loop time-optimal results. The purpose of this section is to devise 
such a strategy. 

The following strategy is based on the piecewis e-line ar model. 
Although further refinement of this model appears necessary, the gen- 
eral form will likely remain intact, and the following will still apply. 

It has been previously noted that the turbine inlet temperature 
limit is an active constraint for the duration of each of the minimum 
time trajectories presented. This fact forms the basis for a closed- 
loop transition control law for fuel flow. Since turbine inlet tem- 
perature is modeled as an output , it can be expressed as a linear com- 
bination of state and control variables, i.e. , 

TIT = c”^x + d^u (6.16) 

For TIT = '^^ax* ^'I’^^.tion (6.16) may be solved for fuel flow (or any 
other control variable for which the coefficient d is nonzero) to 
yield.: 

w = ^ (t _ - c\ - d„A - d„IGV - d.HVS) (6.17) 
f \ max 2 noz 3 ^ / v * 

It still remains to find closed-loop control laws for A , IGV, and 

noz’ * 

HVS. We will consider the case for which IGV and HVS are limited to 
scheduled values ±7 degrees. Examination of a number of optimal tra- 


jectories including those shown in figures 6.3 and 6.4 reveals the 
following facts ; 


(1) always starts at the higher value of 3.2 feet squared, 

then switches to the lower value of 2c8 feet squared. Furthermore, 
the switch occurs at fairly constant values of > averaging about 

11 000 rpm. 

(2) IGV always starts at the lower value (scheduled -7 degrees), 

then switches to scheduled +7 degrees. This switch also occurs at 
fairly constant values of averaging about 12 000 rpm. 

(3) HVS always has the scheduled -7 degrees value, never switching 
to the larger value o 

Based on these observations, the following suboptimal closed-loop 
strategy is suggested. 


noz 


( 


3.2 ft , N £ 11 000 rpm 

’ comp 

2 . 8 ft^ , N £ 11 000 rpm 

comp 


IGV = 


(^Scheduled -7 deg, 
V,.Scheduled +7 deg, 


N £ 12 000 rpm 

comp 

N >.12 000 rpm 

comp 


(6.18) 


HVS = Scheduled -7 deg 


w- = ^ (t - g'^x - d„A - d„IGV - d.HVS) 
f 6.^ \ max 2 noz 3 4 / 

The closed-loop control strategy shown in equation (6,18) was applied 
to the piecewis e-line ar model, and accelerations were obtained for 
various initial values of PLA to intermediate thrust. Results are 
presented in figure 6.13 for an initial PLA of 24 degrees. Terminal 
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error is presented as a function of acceleration time. Similar re- 
sults for optimal open loop accelerations are repeated from figure 6.5 
for comparison „ It can be seen that the optimal and suboptimal results 
are virtually indistinguishable. 

It is also of interest to determine the acceleration time required 
to reduce the terminal error to a fixed value, say 0.05. From fig- 
ure 6,13, it can be seen that the required acceleration time for a 
terminal error of 0.05 is 0.80 second for an initial PLA of 24 degrees, 
for both the optimal and suboptimal results. Similar results were ob- 
tained for other initial values of PLA, and the results are shown in 
figure 6.14c As in figure 6.13, optimal and suboptimal results are 
virtually indistinguishable. 


ACCElfRATION TIME. sec 



INITIAL POWER LEVER ANGLE, PLA, deg 
Figure 6 . 14 . ^ Time required to acceierate to intermediate thrust with a fixed 

terminai error of 0. 05. IGV and HVS, scheduled ±7 deg. 


CHAPTER VII. CONCLUSIONS 


Minimum time accelerations of aircraft turbo fan engines are pre-* 
sented. The calculation of these accelerations is made by using a 
piecewise-linear engine model, and a new algorithm based on nonlinear 
programming. Use of this model and algorithm allows such trajectories 
to be readily calculated on a digital computer with a minimum expendi- 
ture of computer time. 

The new algorithm may be used for solution of optimal control 
problems which are nonlinear in the state variables, and linear in 
the control variables o It is shown that the optimal control for such 
problems is bang-bang, except for possible singular arcs, which are 
not considered. The algorithm requires that a nominal bang-bang solu- 
tion be found that satisfies the system equations and terminal con- 
straints. The Euler-L a grange equations are not utilized in the deter- 
mination of this feasible solution; it is generally not a local mini- 
mumo Equations are derived for the determination of the Lagrange mul- 
tipliers (sensitivity functions) which correspond to the initial fea- 
sible solution. These sensitivity functions are then utilized, along 
with nonlinear optimization (gradient search) techniques , to improve 
the feasible solution. 

The new algorithm has several advantages over methods currently 
in use for solution of such problems. First, the system dynamic equa- 
tions are uncoupled from the Euler-Lagrange equations: the Euler- 

Lagrange equations are not utilized in the determination of the initial 
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feasible solution. Second, use of the new algorithm minimizes the 
number of variables involved in the gradient search. With the new 
algorithm, the search variables are the constraint switching times - 
the times at which switches take place from one set of constraint 
functions to another. Other methods currently in use generally dis- 
cretize the trajectory into a large number of intervals, and search 
for the optimal values of the controls for each interval. 

The algorithm makes use of the fact that the optimal values oi 
the control variables are determined by an equal number of active 
constraints. This situation always exists (except for possible singu- 
lar arcs) when the performance criterion, system equations and con- 
straint equations are all linear in the control variables. However, 
the control may also be fully determined by active constraints for 
some problems for which the equations are nonlinear in the control var- 
iables . It should be possible to extend the theory and algorithm pre- 
sented herein to include such problems; this is an area worthy of fur- 
ther resear’ch. 

The algorithm presented herein is used to find minimum-time accel- 
eration histories for the ElOO engine. A piecewise -linear engine 
model is used, having three state variables and four control variables. 
Minimum time solutions are obtained, and the resulting control his- 
tories are used as inputs to a nonlinear simulation of the FlOO engine 
to verify the accuracy of the piecewise-linear solutions ^ 

Comparison of the nonlinear and piecewise-linear solutions re- 
vealed significant differences in a few of the engine responses, the 
worst being a 35 percent error in the turbine inlet temperature, with 
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several of the variables matching to within 5 percent. To determine 
the source of these errors, the transient response to a small- 
perturbation fuel flow step was calculated using the nonlinear and 
three-and-sixteen-state linear models. It was discovered that there 
were significant differences in several of the engine responses for 
this case as well, notably fan speed and surge margins, for both the 
three- and sixteen-state variable linear models. 

Based on these results, it appears that the calculation of linear 
models from a nonlinear simulation is a difficult task, and that there 
may be significant nonlinearities in the control effectiveness param- 
eters. The linear models used herein were obtained by P&WA by using 
the offset-derivative method. There are at least three methods for 
identification of low-order linear models for aircraft engines which 
merit further study: (1) identification of high -order linear model 

via offset derivative method, followed by modal reduction to low order 
model; (2) least-squares identification of high-order model, followed 
by modal reduction to low order model; (3) least-squares identifica- 
tion of low order model. 

Once accurate linear models are obtained, it is also of interest 
to determine the affect of state and control variable nonlinearities 
on the accuracy of the piecewise-linear results. Further research 
into the identification of simple nonlinear models may also be indi- 
cated. 

Results presented herein indicate that improved steady-state and 
transient performance may be obtained by using optimal control strat- 
egy. Such strategy sometimes calls for operation of the controls in 


124 


a manner which has not been previously tested experimentally, even 
though the nonlinear simulation contains equations for predicting the 
effect of such control action. Further experimental testing is indi- 
cated in order to systematically explore the steady-state and tran— 
effects of of f— scheduled control action, and to determine if 
the predicted performance gains can be achieved. 

In this report, it has been assumed that the control variables, 
i.e. , fuel flow, nozzle area, inlet guide vane position, and compressor 
variable vane position, could jump Instantaneously from one value to 
another. Actually, there are rate limits which apply to each of these 
control variables. It is possible to calculate optimal trajectories 
for which the control rate limits, in addition to the other con- 
straints , are adhered to . This can be done by elevatxng the control 
Yaxiables to the role of state variables, and using the control vari 
able rates as the new control variables. Furthermore, the algorithm 
presented herein can be used to solve this modified problem. 

In addition to the open-loop optimal control strategy derived 
herein, a suboptimal closed— loop control logic is also derived which 
gives close approximation of the open loop performance. However, con- 
clusive evidence on its general applicability would require extensive 
simulation and testing under many different flight conditions . Fur- 
thermore , implementation of the closed-loop control depends ort being 
able to accurately and rapidly sense all engine states, and model all 
outputs. Further analytical research could help to identify types 
and accuracy of sensors needed to accomplish the control objectives. 


appendix a 


OPTIMIZATION WITH STATE VARIABLE INEQUALITY CONSTRAINTS 
In chapter II, the optimal control of a dynamic system which is 
nonlinear in the state variables and linear in the control variables 
is considered, and necessary conditions for optimality are derived. 

In chapter III, a new algorithm for solution of such problems is pre- 
sented. In both chapters II and III, the path constraints are assumed 
to be either control or combined state/control inequality constraints, 
but not state variable inequality constraints. In this appendix, the 
results presented in chapters II and III are extended to apply to 
problems with state variable inequality constraints. To illustrate 
the use of the numerical technique, a state variable inequality con- 
straint is added to one of the example problems from chapter V, and a 

local minimum solution is obtained o 

A path constraint in which the control does not appear explicitly, 
i.e., c.(x,t) £0, is called a state variable inequality constraint. 
When a state variable inequality constraint is active for a finite 
time period (i.e», c^(x,t) =0), the control may be determined from 
the requirement that all time derivatives of c^(x,t) must also be 
identically zero during this time period. If s time derivatives of 

c,(x,t) must be taken before the control appears explicitly, c^(x,t) 

^ _ th 

is called a state variable inequality constraint of order s. The s 

derivative of c^(x,t), i.e., d®c^(x, t) /dt® , then serves as the path 

constraint from which a component of the control is determined. There 
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are s additional equality constraints. 


c^(x,t) = 


dc^(x,t) 

dt 


d® ^c^(x,t) 


dt' 


s-1 


= 0 


which must be satisfied at the initial point of the boundary segment 
on which the control is active. These are sometimes called "tangency'' 
constraints o 

1. Feasible Solution 

In chapter III, a feasible solution having at least p segments 
is obtained for a problem in which there are p terminal constraints. 
If one of the active constraints for a given segment is a state vari- 


able inequality constraint of order s , then there are s additional 
point constraints which must be satisfied at the start of that segment. 
In this event, there must be s additional trajectory segments, whose 
durations provide the additional degrees of freedom necessary to 


satisfy the s constraints. 

2„ Calculation of Lagrange Multipliers 
It is shown in chapter III that once a feasible solution is ob- 
tained, the Lagrange multipliers can be uniquely calculated. This is 
also the case if one of the active constraints is a state variable in- 
equality constraint. Suppose, for example, there is a state variable 
inequality constraint of order s active in phase k. Then there are s 
additional point constraints which must be satisfied at t = 
part of the determination of the feasible trajectory. There are also 
s additional trajectory segments for this case. 

It can be shown (e.g., ref. 13) that at such a point constraint, 
there is a jump in 1 which is proportional to the gradient of the 
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constraint. In this case, we have 


j=l 


J_ 

9x 


-dt" 


T 


(Al) 


For simplicity, we will consider the case in which a single first 
order state variable inequality constraint is active in segment k. 
The results which will be derived can be extended to the case of mul- 
tiple , higher order state variable inequality constraints. For a 
first-order constraint, we have 




"k-1- " ^"k-l^ ' ^ 

Because of the point constraint which must be satisifed at t = 


‘Tc-I’ 

there are (p + 1) trajectory segments, and H must be continuous at 
the p interior points, t t „ , t . In addition, for final 

X p 

time free, equation (3.4) must also be satisfied. 

Except at t = the requirement for continuity of H is 

achieved by satisfying (3.3). However, the requirement for continuity 


of H at t = t. 


must be given special consideration because of 


the jump in A at t = The change in H at tj^ ^ is given by 

AH = + SU(4_J^)) - 




t 


[f + gu(t^_^)J - [f + gu(t^_pi 

+ - u(t- ) ) 






+ ^ I? 
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But for t = t 


k-1^ 


dc 3c r ^ . / + \ 1 


Therefore, we have 


AH = - u(t-.pl at t . 

(A3) 

In order to solve for e and the p values of v , we find 

0 

(p + 2) backward solutions of the A equation (2.13). The first 
(p + 1) of these solutions are identical to (3.5); A^^^^^(t) is ob- 
tained by integrating (2.13) backward, starting at t , with initial 

iC*“X 

conditions A^^ ^\-l^ ~ ” lx *^^k-l^ * superposition, A(t) is 

given by 

A(t) = A^^^(t) + A(t)v + eA^P‘‘'^^(t) (A4) 

The (p + 1) equations for the determination of e and v are given by 

[u(t't') - u(t:)]^[b(t ) + g^(t.)A^°>(t.) + g^(t.)A(t.)v 

+ g'^(t^)A^P'^^\t)t] » 0 

i ~ 1 , • • • j (k — 2) , k , . . . , p (A5) 

+ g-"(t^.^)A(t^_pv - = 0 (A6) 


i(tp - b’’^(t,)D'’^(tpc(tp + Si. (t^) + 




V = 0 (A7) 


The calculation of v and e proceeds exactly as in (3.10) to (3.12), 
and will not be repeated here. 


I 


I 
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3. A Numerical Example 

To illustrate the use of the equations which have been presented, 
a numerical example vAiich includes an active state variable inequality 
constraint will be solved. Consider the problem of finding u(t) 
which transfers the system 

X = -X + u + 0.25 
y = -4y + u 

from initial conditions (xq, y^) = (1, 1) to terminal conditions 
(Xj» Yf) = (0» 0) minimum time. The system is subject to in- 
equality constraints 

u] 1 2 
y + 0 . 2 1 0 

This problem is identical to the one-piece linear approximation 
studied in chapter V, but with the addition of a state variable in- 
equality constraint, y + 0.2 2. 0. 

When the state variable inequality constraint is active, i.e. , 
y + 0.2 = 0, the control is determined from 

y = -4y + u = 0 

and since y = -0.2, we have u = -0.8 as the equivalent path con- 
straint. 

Inspection of the optimal, solution in figure 5.2 shows that the 
constraint y + 0.2 ^0 is violated. Therefore, we look for a fea- 
sible solution in which the constraint y + 0.2 ^ 0 is active on one 
of the segments o Specifically, we assume the optimal control history 
is given by 


(A8) 


(A9) 
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u = “2, 0 1 t 1 

u = -0.8, j 1 t i (AlO) 


u = 2, t 2 £ t 1 

and we seek values of t^, t 2 » and such that the terminal con- 
straints = 0 and the intermediate constraint y(tj^) = -0.2 

are satisfied. Solution of the multipoint boundary value problem 
results in t^ = 0.402, t 2 = 1.002, t^ = 1.086. 

The Lagrange multiplier time history which corresponds to this 
trajectory may be calculated once the three parameters V 2 , and 

e are known. These parameters are determined by requiring that equa- 
tion (A5) (for i = 2), (A6) (for k = 2) , and (A7) be satisfied. The 
resulting equations are given by 

2.25 + 2^2 = -1 

Q.919 +0.715 V 2 = 0 

0.615 + 0.143 - e = 0 

Simultaneous solution of these three equations results in 

- 3.12, V 2 = e = 1.35 

The final conditions on A are therefore given by 

A^(tp = 3.12, Ay(tp = -4.01 

Also, A is discontinuous at t, : 

’ V 1 


Xy(tp = Ay(tp + 1.35 

The switching function (X + X ) is presented as a function of time 

X y 

in figure Al. It confirms the fact that the initial feasible solu- 
tion is a local minimum. The optimal trajectory is shown in figure A2. 
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Fhnire A-1 - Switching function and control profile lor optimization 
nrohlem With State variable inequality constraint. 
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appendix b 

derivation of LAGRANGE MULTIPLIER FUNCTIONS FOR lEASIBLE SOLUTION 

Necessary conditions for the optimal control of a dynamic system 
which is nonlinear in the state and linear in the control are derived 
xn chapter II. It is shown that except for singular arcs, the optimal 
control (r variables) is always determined by r active constraint 
boundafi^s. In chapter III, the nature of the optimal control strategy 

xs used as the basis for a new algorithm for the solution of such opti- 
mization problems . 

In chapter III, the first step is to obtain a feasible solution 
which satisfies all path and terminal constraints. This is accom- 
plished by dividing the trajectory into at least p segments (for a 
problem with p terminal constraints) and requiring the control to 
be determined by a set of r active constraints for each segment; p 
of the junction times are varied in order to satisfy the terminal con- 
straints, The Euler-Lagrange equations are not utilized in the deter- 
mination of this feasible solution; it may or may not be a local mini- 
mum. 

Once such a feasible solution is obtained, it is useful to cal- 
culate the Lagrange multiplier time histories, functions of which are 
switching functions which predict the change in the performance index 
obtaxnable from small changes in the control histories. In chap- 
ter III, xt xs shown that these Lagrange multipliers can be easily 
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, , d The calculation scheme makes use ol the di£- 

and uniquely oalcu a e . continuity of the 

• for the multipliers, and the 

ferential equatxo ntiolier differential 

. . function points. The Lagrange multiplier 

Hamiltonian at 3n caquirement for con- 

. are derived in this appendix. Also, the 

, section points is proved here . 
tlnulty of the Hamiltonian a 1 ^ 

cedure followed in chapter 

d lh these segents are used in conjunction 

suitching times assoc 

a -set search technique to systematically 
with a gradient nartial derivatives 

. dex The gradient search technique requires partia 
ance index. switching times; 

of the performance index with respec 

^or-Tved in this appendix, 
r-it,! derivatives are derived m 

these partial problem is as 

follows. It is desired to find the values of 
Uhich result in minimising the performance index 


J (t)(x^,tp + 


(a + b'^u)dt 


subject to the system dynamical equations 


X = £ + 


and terminal constraints 

if.(Xj,tf) - »• ’ 

to addition, the control must satisfy 

u = -b:\, '^1-1-' -"d 


i P 


= 1, 2, . . “ ip 




I 


J 


1. 
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w 


•\ 



J" 


i=l 


Iraa . Sb’^ + jT /3£ ^ \ I / °°1 - I J-I I ;il 

S |_sr toT “ ■^ '■i “j ■ '=1 \ 3x + “l 3x ) + 


6x 


+ (X^ + pj + b^)6u dt + j(a + b^u)dt - xT 6:^^^ 

i=l 

Since u is constrained by (B4) , the 6u are not free but are 
given by 

(^^7 -T 


(B7) 


(B8) 


If (B4) and (B8) are substituted into (B7) , and we choose ^^(b) such 
that 


IN 

X 


T 9C, 

|a^|^^:Te. i C. +bV"^ 

9x 3x X 1 3x X X dx 


/ 9dT^ 

i V 9x 9x i i ® 9x i 


- gD 


rr, 3C.^ 

-T X 

i 9x 


(B9) 


then (B7) simplifies to 


^ ^ i=l 


T . 


[(a + b^u)dt - X^ 6x]^-" 


(BIO) 


i-1 


Note that (B9) can be written 


i 


k 


I 
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(a - bD-\) - Xj (f - gDf C,) 


(BID 


which Is identical to (2.13). 

The variations 6x which appear in (BIO) apply for fixed time 

t and can be expressed in terms of the total variation dx and 
i’ ■ 

variation of time dt. 

fix = dx - X dt (B12) 

Now if we use (B12) in (BIO) , and choose 


f 3xg 3x^ 


(B13) 


and if t^ is free. 


H(tD A a + b"^u + x'^(f + gu) = 


34> „T H 

3t. 3t. 

f f 


(B14) 


then (BIO) simplifies to 
w-1 


6J 


■E 


(a + b”^u + /\Tx)dt - xT dx 

X X 


■^t. 

X 


i=l 

w-1 


t. . 
x-1 


w-1 


i=l ” i=l 


b'^(t.)u(t^) + xj(t^x(t") 


- b'^(t,)u(tt) ~ xL (t,)x(t^)^dt. 


X- -'D '^i+l'*x'‘‘'"^Dl i 

Note that (B13) and (B14) are identical to (2,7) and (2.8). ¥e can 
also choose 


(B15) 


so that the variation of J becomes 




I 
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6J 




i=l 

w-1 


.y^[H(tp - Hc4)]dt, 
i=l 


(B16) 


If a particular t^ is one of the free junction times, the coeffx 


cient of dt^ mus 


t vanish for a local minimum solution. 




H(t7) = H(t ) 
1 1 


(B17) 


for free junction times t^. On the other hand, xf t^ t^^ is one 
of the fixed times, then we have 


= ii— = H(t“ ) - H(t^ ) 
3t 9t sw sw 

sw sw 


(B18) 


APPENDIX C 

MODEL ORDER REDUCTION 
Suppose we have a linear system 

X = Fx + Gu 

subject to linear inequality constraints 

Cx + Du + E 1 0 

Also, suppose there are n states, and n^ of the eigenvalues of F 
are much smaller than the remaining n^ = n - n^. Then the dynamics 
associated with the larger-eigenvalue modes are much faster than 

the dynamics of the n^^ lower eigenvalue modes. Therefore, the 
higher frequency modes will nearly always be in equilibrium (or oscil- 
lating at high frequency about equilibrium). In such a case, the cal- 
culation of optimal trajectories can be made significantly easier if 
the high frequency modes are assumed to be always in equilibrium. In 
carrying out that approximation, it is also desirable that the low fre- 
quency eigenvalues be retained exactly. Furthermore, it is desirable 
that n^^ of the original states be retained in the lower order approx- 
imate model. 

Specifically, we seek a low frequency approximation of the system 
in which the first n^ states of the original system are retained, 
and the eigenvalues of the reduced order system are exactly equal to 
the n^ smaller eigenvalues of the original system. We proceed as 
follows. 
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V 


We can find a transformation T which block-diagonalizes the 
system: 

X = Tm; t"^x = m (Cl) 

where the m are modal coordinates. Differentiating the above gives 


„-l. 


m = T "^x = Am + Su (C2) 

where Z t and A = T"^FT is block diagonal. If (C2) is par- 

( 1 ) 

titioned into low frequency modes m and high frequency modes 


m^^^ , the result is 




The quasi-steady approximation is 

.( 2 ) . „ 

m =0 


(C3) 


which results in 

„(2) . (cA) 

/T \ 

and denote the 

states of most interest in the low frequency approximation. We have 


-1 


Now let x^^^ have the same dimension as m^^^ 


( 1 ) _ ( 1 ) . ( 2 ) 


Differentiation results in 




( 1 ) 


= T 

11^^ 11 


" ^12"^ 1 




^ ll^'- 11 

x^^^ H 

y Ft 

il'’ 


+ 

-1 ( 2 )“^ ( 2 ) 
■^^llWl2^ 


u 


fns') 


/ 
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APPENDIX D 

THREE-STATE LINEAR MODELS. SEA-LEVEL-STATIC 
(a) PLA = 20° 


-2,09 

.0944 

U.6 

''7.27 

.756 

V7.38 


/ 0.239 
/ 1.77 
' .0521 

- 2.01 
i-16.1 
\ .825E-2 

\ .0233 

/ 0.314 
.445 
,436 
12.3 
222 

y-0,621 

\6.09 


2.09 

-0.251 

-0.457 

0.054 

-0.0105 

-1.87 


-0.416E-2 

.0839 

-0.587 

3.51 

37,9 

-1.60 

2.06 

-0.270E-2 

.515E-2 

-0,0130 

.150 

1.98 

-0,0155 

.0382 


-0.537 

-0.163 

-18.9 


-0.478E-2 

.363E-3 

,0584 


0.336 

-0.108 

.0245 

-3.55 

-5.24 

.384 

•189E-2 

0.150E-2 

.0114 

-0.508E-3 

-0.376 

-0,0285 

-0.199E-2 

.108E-2 


38.79 

91.92 

161 


-0.00398 

-0.00109 

■0.00378 


11.34 

300 

-250 

-3050 y 


11.26 
73.3 
54.1 
102 
83.8 
120 
^ 73.2 


-0 . 787E-4 
-0.131E-3 
.237E-2 
-0.476E-2 
-0.491 
.686E-2 
-0.463E-2 
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I 
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(b) PLA 


= 36' 



0.905 


0.988 

-0.232 

o435 

0.0138 

-0.0133 

-4.67 


-0.753E-3 


-0.462\ 
-0.106 
-7.91 / 

-0.0350 

-0.00478 

.370 


0.243 


/ 77.2' 
X = ( 109 

® \226 j 

-0.0021l\ 
-0.00139 
.0135 / 


/ 52.2 


/ 3.10 

.0841 

-0.107 \ 

/ 

1 -0.102 

-0.135 

-0.0137 

/ 

4.39 

1.73 

-1.89 


\ -5.63 

28.0 

-0.594 

\ 

\-0.177 

-0.398 

-0.00598y 

\ 

\l.05 

1.14 

.137 / 



u = 



0.285 

-0.0175 

0.0130 

0.117E-4 

.143 

.108E-2 

.0439 

-0.327E-3 

.703 

.291E-2 

-0.193E-2 

.102E-2 

2.27 

.0589 

-0.199 

-0.541E-2 

1.27 

-0.0720 

-0,0328 

-0.118 

1.46 

.0123 

-0.383E-2 

.261E-2 

2,05 

-0.862E-2 

.0136 

-0.370E-2 
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(c) PLA = 52° 


( 


-2.27 

-0.498 

22.8 


( 


2.00 

1.30 

7.41 


2.23 -0.520\ 

-1.96 -0.0639) 

.936 - 6.11 / 



-0.0366 -0.0434 

.00653 -0.0103 

-6.24 .541 


-0.531E-3 \ 
-0.445E-2 ) 
-0.166E-2 / 


u 

e 




0.0466 

.0937 

-0.136 

2.59 

22,3 

-0.321 

2.21 


0.224 \ 
-0.074l\ 
-0.0249 \ 
- 1,86 
-0.399 
- 0.0596 
1.91 / 



/ 


3.34 

-0.0147 

0.0240 

0.675E-4 

.256 

.0163 

.0665 

-0.550E-3 

.630 

.977E-4 

-0.528E-2 

.591E-3 

2.82 

.0244 

-0.190E-2 

-0.634E-2 

1.29 

-0.115 

-0.765E-2 

-0.0906 

1.36 

-0.0107 

-0.0102 

.162E-2 

^ .962 

-0.268E-2 

.03C7 

-0.138E-2 
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(d) PLA = 67° 


F = 


-3.37 

-0.130 

-21.3 


3.13 

-2,87 

1.93 


-0.493\ / 94. 4\ 

-0.072 ) X = I 122 ) 

-6.03 / ® \330 / 


/ 1,59 -0.0354 

G = I 1.17 .0165 

\5.74 -7.34 


-0.0848 0.155E-2\ 

-0.795E-2 -0.574E-2 ) 

.804 .0173 / 


u 


e 



C = 



0.0511 

.0249 

-0.187 

2.00 

12.3 

-0.465 

3.16 


0.219 \ 
-0.0198\ 
-0.0353 \ 
-1.42 
-0.245 
-0.0817 
.282 I 



D = 



-0.00523 

,0432 

-0.294E-2 

.135 

.206 

-0.025 

,0335 


0.0342 

.0799 

-0.191E-2 

-0.671 

-0.0967 

-0.437E-3 

.0215 


0.807E-4 
-0.561E-3 
0.436E-3 
-0 .444E-2 
- 0.110 
0.119E-2 
.819E-3, 
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APPENDIX E 


COMPUTER PROGRAM LISTINGS 


MAIN PROGRAM 
IMPLICIT REAL*8 


(A - H, 0 - Z) 


COMMON/CA/C (8,3, 5) ,D (8,4,5) ,E (8) ,XI (3,5) ,01 (4,5) ,TI (8,5) ,X0 (7) 


COMHON/CB/A (7,7,20) ,BX (7 ,20) ,T (20) , VAL(80) ,JH (80) , HIT(20) ,JHIT(20) 
COMMON/CO/ALFA (7,20) , BETA (20) , FARMS (3,4) ,TS(2 0) ,XSTO(7,20) , 

1XLSTO (7,20) ,XLO (7) ,IT (3) , JM (5) ,NPHASE 
COMMON/MTRANS/USLOPE (2) ,NN 
COHMON/NNOUT/NOUT 
COHMON/NS/NSTATE 
COMHON/OUTCNT/NOUTLM 

COMHON/TTRANS/PAR (20,20) ,PARIN(3,3) ,PAES(3,3) ,X(7) ,KFLAG,KNAX, 
1KOPT(20) ,KOUNT,NOPT (20) ,NOPTS 
DIMENSION Z(5) , ZL(10) ,ZOPT(5) 

EXTERNAL FUNGRD 
LOGICAL CONV,UNITL, LOCAL 


NSTATE=7 

IF(NOPTS.EQ.O) GO TO 11 
DO 10 I=1,NOPTS 
J=NOPT (I) 

10 ZOPT(X)=TS (J) 

11 CONTINUE 


NOUTLH=-1 
IPRINT=-1 
12 CONTINUE 

NLDIM=NOPTS* (NOPTS- 1 ) /2 
IF (NLDIM.EQ.O) NLDIM=1 
KOUNT=0 


NN=0 
NOUT=0 
N=NOPTS 
DO 1 1=1,20 

1 KOPT(I)=I 
CALL SETUP 

zx=o. 

ZOL=1 . D-5 
ZPSMCH=16.0»* (-13) 

ZTA=0.5 

ZTEPMX=1.0 

UNITL=-TRUE. 

LOG AL=. FALSE. 

IF (NOPTS.EQ. 0) GO TO 3 

CALL QNHDEH (N, NLDIM,ZOPT,ZX , FUNGRD, ZL, Z,ZOL,ZPSMCH,ZTA, NFTOTL, 

1 NGTOTL, NITER, ZTEPMX,0 NITL, LOCAL, IPRI NT, CON V) 

CALL OUTPT 
GO TO 2 

3 T(1)=TS(1) 

DO 4 I=2,NPHASE 

4 T(I)=TS(I)-TS (1-1) 

CALL TRAJ 

CALL LCALC 
CALL OUTPT 

2 CONTINUE 
STOP 
END 


I 


li|8 


11 

3 

1 


14 


12 


21 


SUBROUTINE FUNGPD(M,ZOPT, IPLAG, ZP,Z) 

IMPLICIT REAL*8 (A - H, 0 - Z) 

R5AL*8 IS 

»E(8) , XI (3,5) ,01(4,5) ,11(8,5) ,X0(7), 

COBHON/CO/ALFA (7,20) , BETA (20) , PARRS (3,4) ,TS (201 . XSTO <7 20 \ 
1XLSTO (7,20) ,XLO{7) ,IT(3) , JR (5) , NPHASE rISTO(7,20) , 

COHHON/LTRANS/ZLAH (7) ,IS (7,7,20) 

COHMON/NNPOT/NOUT 

COMHON/NS/NSTATE 

CORBON/TTRANS/PAR(20,20) ,PARIN(3,3) , PARS (3,3) , X (7l .KPLAG KB AX 
1KOPT(20) ,KOUNT,NOPT(20),NOPTS r A (/), KPLAG, KRAI, 

DATA^XER/7*o” (7) ,TERP(7) ,Y (7) , Z (N) , ZOPT (N) 

DO 11 I=1,NOPTS 
J=»OPT(I) 

TS(J)=ZOPT(I) 

T(1)=TS(1) 

DO 1 1=2, NPHASE 
T(I)=TS(I)-TS(I-1) 

CALL TRAJ 
T(1)=TS(1) 

DO 5 1=2, NPHASE 
T(I)=TS(I)-TS(I-1) 

KRESET=0 

IF(KRESET.EQ.I) KRESET=2 
DO 14 1=1, NPHASE 
IP(T(I).GE.0.0)GO TO 14 
CALL RESET (I) 

KHESET=1 

CONTINUE 

IF(KRESET.EQ. 1) GO TO 7 
IP (PRESET. EQ. 2) CALL TRAJ 
IP(IPLAG.EQ. 2) GO TO 12 
ZP=0.0 
DO 2 1=1,3 

ZF=ZP> ( (X (I) -XD (I) ) /XD (I) ) ♦♦2/2. 0 
IP(IPLAG.EQ. 1) RETURN 
CALL LCALC 
JJ=1 

KX=K0PT(1) 

LX=NOPT(KX) 

RX=K0PT (NOPTS) 

LOPT=NOPT (MX) 

DO 20 I=1,LOPT 
IF (I. NE. LX) GO TO 20 
DO 21 J=1,NSTATE 
Y(J)=BX(J,LX)-BX(J,LXf1) 

DO 21 K=1,NSTATE 

IX-I+ 1 ~ K,LX*1) -A (J , K, LX) ) ♦XSTO (K, I) 

IF ( IX. GT. NPHASE) GO TO 2 5 
DO 22 J=IX, NPHASE 
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23 


24 
22 

25 


20 


« 1 XES.HST ATE# 1) 

CALL WPDATE(Y,T(J) .^^^^^2'^ ' 

IF (JHIT (J) -EQ-0)^0 * 

DD=0.0 
DN=0.0 

DO 23 /R jY 

s5:s5:r«trK'/;,'.“^s^ ’ 

rD.S»»«UK!”*MK,I..3)*XST0aW) 

PAEIJ,!)--”"/”” 

(K , a. 1 ) -B« ( K. j) I *» <•> ' 

r,Kr-TrK;-"(MK!i^*i. -MK.i..^) > 

CONTINUE 

KX=K0FT(JJ) 

LX=N0PT (KX) 

CONTINUE 

JU=1 

KX=K0PT (1) 

LX=N0PT (KX) 

DO 4 I=1#LOPT 
IF (!• EE. LX) GO TO 4 
2(KX)=0-0 

IX=I+1 

IF (IRON. EQ. 1) GO “^0 
IF ( IX . GT . NP H ASE) GO 
DO 28 J^IXpN PHASE 

IF(JHIT (0) .n#-iTV 

Z(KX)-Z(BX);“ST0( 

27 CONTINUE 

28 CONTINUE 

rEHM6!lB)B<>‘'i(XX) ,TB,Z(KX) .BX 

aa=u3+'’ 

kx=kopt(JO) 

LX=N0PT(KX) 

CONTINUE ir7(\ 91 I2) 

FORHAT (1H »3(G20.9)#I2) 

IF(NOUT-EQ.O)CALL ootp 

UO0T=1 

RETURN 
END 


8 


26 

TO 


26 


TO 28 


4 

10 


I 


V 
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SUBROUTINE LCALC 

IMPLICIT BEAL*8 (A - H, O - Z) 

EEAL*8 ID,rS,ISIN 

^COMMON/CA/C(8,3,5) ,D (8,4,5) ,E(8),XT(3,5) , UI (4 , 5 ) ,TI (8, 5 ) ,X0(7), 

COWMON/CB/A (7,7, 20) ,BX(7,20) ,T(20) , VAL(80) , JN (80) , HIT (20) , JHIT(20) 
COMMON/CO/ALFA (7 ,20) ,BETA(20) , FARMS (3,4) ,TS (20) ,XSTO (7,20) , 
1ZLSTO(7,20) ,ZL0 (7) ,IT(3) ,JM (5) ,NPHASE 
COMMON/CX/F (3,3,5) ,G(3,4,5) ,DD(4,4,20) , ID (7,7) 

COMMON/EPSS/EPS (20) 

COHMON/LTKANS/ZLAH (7) ,IS (7,7,20) 

COMHON/NS/NSTATE 

DIMENSION GSP(7) ,GSM(7) ,TEMP(7) ,ZER (7) ,ISIN (7,7) 

DATA ZER/7*0.0/ 

DO 2 1=1,3 

ZLAM(I)= (XSTO(I,NPHASE)-XD(I) )/XD(I) **2 

2 ZLSTO (I,NPHASE)=ZLAM(I) 

DO 11 I=4,NSTATE 
ZLAM(I)=0.0 

11 ZLSTO (I,NPHASE) =ZLAM (1) 

DO 3 I=1,NPHASE 

IA=NPHASE-I+1 

TX=-T(IA) 

CALL UPDATE(ZLAM,TX,A(1,1,IA) ,ZER,NSTATE, 2) 

IF(IA.EQ.1)GO TO 3 
DO 4 J=1,NSTATE 

4 ZLSTO (J,IA-1)=ZLAM (J) 

IF(JHIT (IA-1) .EQ.O) GO TO 3 
DO 5 J=1,NSTATE 
GSP(J)=BX(J,IA) 

GSM (J)=BX(J, IA-1) 

DO 5 K=1,NSTATE 

GSP (J) =GSP(J) <-A (J,K,IA) ♦XSTO(K, IA-1) 

5 GSM(J)=GSM (J) +A (J,K,IA-1)*XSTO (K,IA-1) 

DOT=0.0 

DOTIN=0-0 

DO 6 J=1,NSTATE 

DOTIN=DOTIN + ALFA(J,IA-1) *GSM(J) 

6 DOT=DOT+ZLAH (J) ♦ (GSM (J) -GSP (J) ) 

EPS (IA-1) =DOT/DOTIN 

DO 8 J=1,N STATE 

ZLAH(J)=ZLAM (J)-EPS (IA-1 )♦ ALFA (J, IA-1) 

8 ZLSTO (J, IA-1) =ZLAM(J) 

3 CONTINUE 

DO 10 I=1,NSTATE 
10 ZLO (I) =ZLAM(I) 

RETURN 

END 


I 


t 
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SUBROUTIKE OUTPT 

IMPLICIT RERL*8 (A - H, 0 - Z) 

D (8.» ,5, ,E (8, , XI (3,5) ,01 (4,5) ,YI (8,5) ,X0 (7) , 

’co5bO»/CB/A (7,7,20) ,BX H -20) ,T (20) . 

COMMON/CO/ALFA(7,20) ,BETA(20) ,PARMS(3,4) ,TS (20) ,XSTO(7,20) , 

1ZLSTO(7,20) ,ZL0(7) ,IT(3) ,JM{5),NPHASE 

COHMON/GX/F (3, 3, 5) , G (3, 4 , 5) ,DD (4,4, 20) , ID (7,7) 

COMMON/EPSS/EPS (20) 

COMMON/LTRANS/ZLAM (7) ,IS (7,7,20) 

COMMON/MTRANS/USLOPE(2) ,NN 

COHMON/NS/NSTATE 

COMMON/OOTCNT/NOOTLM 

DIMENSION S(5) ,X(7) ,OUT (8) , KAPPA (4) ,TEMP(7) ,ZER(7) 

DATA ZER/7*0./ 

IF (NOUTLM.LT.O) RETURN 
M=1 

IM=JM (1) 

DO 50 1=1,4 
IF (IM.GT.O) GO TO 50 
M=M*1 
IM=JM (M) 

50 CONTINUE 
KK*1 

WRITE(6,7) TS(20) 

WRITE (6,13) 

WRITE (6, 14) 

WRITE (6,20) 

13 format r* I' ',4 x','tiL'',71!*4hRUSTS5X, AIRFLOW', 4X,'TITS0X,'SMPAN', 

16X *SHCOM ' ^ 6 X, ' FTIT ’ ,7X , 'PT3 ' ,8X, ' NFAN * ,7X, 'NCOHP' , 6X, *PT7H ) 
pShMaHL •W^,9X,•ANO^,7X,WG^,8X,•HVS•,8X,•KAPPA1•,5X, 
VxjrpL-'"5L"KAPP;^%X,'i;AP««',5I,'H',10X,'Dl',9X,'D2') 

20 FORMAT (5X , • D3 ' , 9X , • D4 • , 9X , ' D5 • , 9X, * TERHEFF ) 

DT=0.1 

TIHE=0.0 


N=1 

DO 41 I=1,NSTATE 
ZLAM (I) -ZL0 (I) 

41 X(I)=X0(I) 

GO TO 16 ^ . 

8 IF (TIME+DT.LE.TS (N) , I ND. KK- EQ. 1) GO TO 1 

GO TO (10,11,18,12) ,KK 

10 TSTO=TIME*DT 
DTSTO=DT 

60 DT=TS(N)-TIHE 

KK=2 
GO TO 1 

11 N=N*1 

IF (N. GT. NPHASE) RETURN 
IF(JHIT(N-1).EQ.O)GO TO 40 
IF(N.NE. IH*1) GO TO 51 
M=M*1 
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iH=Jn (M) 

51 DO 53 J=1,NSTATE 

53 ZLAM(J)=ZLAM (J) +EPS (N-1 ) ♦ ALFA (J, N-1 ) 
40 DT=0.0 
KK=3 


GO TO 16 
18 DT=TSTO-TIME 

IF(TS (N) .LT.TSTO)GO TO 60 
KK--4 
GO TO 1 
12 KK=1 

DT=DTSTO 
GO TO 8 

1 CALL DPDATE(X,DT,A(1,1,N) ,BX(1,N) ,NSTATE,1) 
TI«E=T1ME*DT 


16 NST=NSTATE-4 
DO 3 1=1,7 
OUT(I)=YI(I,M) 

DO 34 J=1,4 

34 OUT (I) =OUT (I) +D(I, J,M)*(X(J + NST)-UI (J,M) ) 
DO 3 J=1,NST 

3 OUT (T) =OUT (I) +C (I, J,M) * (X(J) -XI (J,M) ) 

TX=DT 

IF(DT.EQ.0.0.OR.TIME.EQ.0,0)GO TO 17 
CALL UPDATE (ZLAM,TX, A (1,1, N) ,ZER, NSTATE, 2) 

17 DO 9 1=1,4 
KAPPA (I) =0.0 
DO 9 J=1,4 
JX=J+NST 

9 KAPPA (I) =KAPPA (I) -DD ( J, I, N) ♦ZLAN (JX) 

H=0, 


DO 15 1=1, NSTATE 
H=H+ZLAM (I)*BX(I,N) 

DO 15 a=1, NSTATE 

15 H=H+ZLAW (II^A (I,a,N) *X (J) 

DO 37 1=1,5 

S(I)=0.0 

DO 37 0=1, NST 

37 S (I) =S (I) + (X (J) -XI ( J,I) ) **2/ (XI (J,1) -XI (J,2) ) ♦♦2 
DO 52 1= 1 , 5 

52 S (I)=DSQET (S (I)/3.0) 

ZF=0.0 


4 

6 

19 


DO 2 1=1 ,NST 

ZF=ZF+((X(I)-XD(I) ) /XD(I) ) *♦2/2.0 
TE=DSQRT (2.*ZF) 

WRITE (6,19)TIME, (OUT (I) ,1=1 ,7) , (X(I) ,1=1,3) , , 

WRITE (6,4) (X (1+3) ,1=1,4) , (KAPPA (I) ,1-1,4) ,H, (S (I) , 1-1,2) 
WRITE (6,6) (S (I) ,1=3,5) ,TE 
FORMAT (3X, 11 (Gil. 4) ) 

FORMAT(3X,3 (G11.4) ,G11. 5) 

FORMAT (3H0 , G1 1. 5, 10 (Gl 1 . 4) ) 

GO TO 8 
END 
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10 


11 

12 


SUBRODTINB RESET (N) 

IHPLICIT REAL*8 (A - H, 0 - Z) 

CONHOM/CB/A (7,7,20) ,BX (7,20) ,T(20) ,?AL(80) (80) ,HIT(20) ,JHIT(20) 

COMNON/CO/ALFA(7,20),BBTA(20),PARHS(3,4),TS(20),xbo(7,20)r 

1XLSTO (7,20) ,XLO(7) ,IT (3) ,J« (5) , NPHASE 

COMMON/TTRANS/PAR(20,20) ,PAHIN (3,3) ,PAHS (3,3) , X (7) ,KPLAG.KHAX, 
1KOPT(20) ,KOUNT,NOPT(20),NOPTS » » » *■ 

DO 12 J=1,4 
H=4*(N-1)fJ 

IF(JN (H) ,NE. JN(M-4) ) GO TO 1 0 
JNT=JN(M+4) 

VALT=VAL (n*U) 

IF(JN (B) .ME. JN(M*4) )GO TO 11 
JMT=JN(B-4) 

VALT=VAL (M-4) 

JN(B)=JNT 
VAL (M) =VALT 
DO 7 1=1,7 
ALFT=ALPA (I,N) 

ALFA(I,N)=ALPA (I,M-1) 

ALFA(I,N-1)=ALFT 

BET=BETA(N) 

BETA (M)=BETA (N-1) 

BETA(M-1)=BET 
H1TT=HIT (N) 

HIT (N)=HIT (N-1) 

HIT (N-1) =HITT 


JHITT=JHIT(N) 

JHIT(N)=JHIT (N-1) 
JHIT(N-1)=JHITT 
TT=TS (N) 

TS (N)=TS (N-1) 

TS (N-1) =TT 
T(N)=TS(N)-TS(N-1) 

T(N+1) =TS (N+1) -TS (N) 
IF(N.GT.2)T(N-1)=TS(N-1) -TS (N-2) 
IF(N.EQ.2)T(1)=TS(1) 

DO 2 1=1,5 
K=JH(I) 

IP(K. EQ. N) JH(I)=K-1 
IF(K.EQ.N-I) JH(I) = K+1 
2 CONTINUE 
DO 3 1=1,20 
K=NOPT(I) 


IP (K. EQ. N) NOPT (I) = K-1 
IF (K. EQ.N-1) NOPT(I) =K+1 
3 CONTINUE 
DO 4 1=1,19 
KX=KOPT (I) 

KY=KOPT(I*1) 

IF (NOPT(KY) .GE.NOPT (KX) ) GO TO 4 
KT=KOPT(I+1) 

KOPT(I*1)=KOPT (I) 

KOPT(I)=KT 



15 ^ 


4 CONTINUE 
CALL SETUP 

¥KITE(6,5)N, (JM(I) ,1-1/5) 

5 FORMAT {1H ,7 (13) ) 

RETURN 

END 
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100 


27 


41 


40 


SUBROUTINE SETUP 

IMPLICIT SEAL*8 (A - H, Q - Z) 

COMMON/CA/C ^ »5) , D (8, 4 , 5) ,E <8) , XI (3 , 5) ,01 (4, 5) ,YI (8,5) ,X0 (7) , 

^COHMON/CB/A (7,7,20) ,BX(7,20) ,T (20) , V AL (80) 

COMMON/CO/ALFA(7,20) ,BETA(20) ,PARHS(3,4) ,TS(20) ,XSTO (7,20) , 
1ZLSTO(7,20) ,ZL0 (7) ,IT(3) ,JM(5),NPHASE 
COMMON/CX/F (3,3,5) ,G (3,4,5) ,DD (4, 4,20) , ID (7,7) 
COHMON/MTRANS/OSLOPE(2) ,NN 

C2(»,a),CS(2,3), ,1(3.51 ,««<*) 

DO 100 1=1,7 
DO 100 J=1,20 
ALFA(I, J)=0.0 
DO 27 1=1,2 
DO 27 J=1,3 
CS(I,J) = 0.0 
CS (1,1)=23./6400. 

CS (2,2)=34./4100. 

DO 41 K=1,20 
DO 41 1=1,4 
DO 41 J=1,4 
DD(I, J,K)=0.0 
NST=NSTATE-4 
DO 40 1=1, NST 

WI (I,1)=1.5*XI (I,1)-0.5*XI (1,2) 

DO 40 J=2,5 
WI (I, J)=xi (I,J) 

M=1 

IM=JM(1) 

DO 50 1=1,4 

IF (IM. GT.O) GO TO 50 

M=M+1 


IH=JH (M) 

50 CONTINUE 

DO 1 K=1,NPHASE 
IF(K.LE. IM)GO TO 2 
M=M + 1 
IH=JM(H) 

2 DO 3 1=1,4 
EE (I) =0- 0 
DO 3 J=1,4 
ci(i,a)=o-o 

C2(I,J)=0.0 

3 DD(I,J,K)=0.0 
DO 4 1=1,4 
LL=4* (K~1)+I 
N=aN(LL) 

NT=0 

IF (N.GT.8) GO TO 5 
19 DO 6 J=1,4 

DD(I,J,K)=D(N,J,M) 
DO 6 L=1,NST 
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6 C2{I,J)=C2(I,J)*C(M,L,H) ♦G(L,J,H) 

DO 7 J=1,NST 

DO 7 L=1,NST 

7 Cl (I, J}=C1 
DO 8 J=1,MST 

8 EE (I) =EE (1)-C1 (1,J) *X1 (J,H) 

DO 9 J= 1,4 

9 EE (I) =EE (I)-C2 (I,J) *UI (J,H) 

GO TO 10 

5 IP (H.GT, 12) GO TO 11 
DD(I,M-8,K) =1,0 
GO TO 10 

11 IP(N.GT. 16) GO TO 12 
DD (1,N-12,K)=1.0 
EE(I) =-TAL(LL) 

GO TO 10 

12 IP(N. GT. 20) GO TO 13 
DD (I,N-16,K)=1 .0 
EE(I)=VAL(LL) 

GO TO 10 

13 IP <N.GT. 22) GO TO 14 
DD (I,N-18,K) =-1.0 
DO 15 J=1,4 

DO 15 L=1,NST 

15 C2{I,J)=C2(I,J) +CS(N-20,L) *G(L, J,M) 
DO 16 J=1,NST 

DO 16 L=1,NST 

16 Cl (I,J) =C1 (I,J)+CS (N-20,L) ♦P(L,J,M) 
DO 17 J=1,NST 

17 EE (I) =EE (I)-C1 (I,J) ♦Xl (J,H) 

DO 18 J=1,4 

18 EE (I) =EE (D-C2 (I,J) *01 (J,H) 

GO TO 10 

14 N=3 
NT=1 

GO TO 19 

10 IF(NT.EQ.I) EE(I)=EE(I) -VAL (LL) 

4 CONTINOE 

CALL DHINV(DD(1,1,K) ,4,DET,LX,HX) 

DO 20 1=1, NST 
DO 21 J=1,NST 

21 A (I,J,K)=P (I, J,M) 

DO 20 J=1,4 
JX=J+NST 

20 A(I,JX,K)=G(I, J,M) 

DO 22 1=1,4 
IX=I+MST 
BX(IX,K)=0.0 
DO 22 J= 1,4 

22 BX(IX,K)=BX(IX,K)-DD (I,J,K) ♦EE(J) 

DO 23 1=1, NST 

BX(I,K)=0.0 
DO 24 J=1,NST 

24 BX(I,K)=BX (I,K)-F(I,J,H) <‘XI (j,H) 
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DO 23 J=1,4 

23 BX(I,K)=BX(I,K)-G(I,J,H) *UI (J,H) 

DO 25 1=1,4 
DO 25 J=1,NST 
IX=I*NST 
A(IX,J,K) =0.0 
DO 25 L=1,4 

25 A (IX, J,K)=A (IX, J,K) -DD(I,t,K) *01 (L, J) 

DO 26 1=1,4 

DO 26 J=1,4 
IX=I*NST 
JX=JfNST 
A(IX,JX,K)=0.0 
DO 26 L=1,4 

26 A(IX,JX,K)=A (IX^.1X,K)-DD(I,L,K) *C2(L,J) 
1 CONTINUE 

B=1 

IK2=JH(1) 

DO 28 1=1,4 
IF(IH.GT.O) GO TO 28 
M=I! + 1 
I«=JH(M) 

28 CONTINUE 

DO 29 I=1,NPHASE 
IF (I.LE.IH)GO TO 31 
M=H + 1 
I11=JM (H) 

31 JH=JHIT(I) 

IF (JH.EQ.O) GO TO 29 
IF(JH.GT.8) go TO 30 
BETA(I)=HIT(I)-YI(JH,f!) 

DO 32 J=1,NST 
ALFA(J,I)=C(JH,J,«) 

32 BETA(I)=BETA (I) +C ( JH, J, H) *XI ( J, M) 

DO 33 J=1,4 

JX=NST*J 

\LFA (JX,I)=D (JH,J,M) 

33 BETA(I)=BETA(I) +D ( JH, J, H) *UI ( J, M) 

GO TO 29 

30 IF(JH.GT.12)GO TO 34 
ALFA (JH-5,1) =1.0 
BETA(I)=HIT(I) 

34 IF (JH.LT.21) GO TO 29 
IF (JH- GT.22) GO TO 3 5 
BETA (I)=HIT (I) 

IF(JH. EQ.22)GO TO 101 

BETA (r)=25. + 23. *4000./6400.-HIT(I) 

GO TO 102 

101 BETA (I)=30.+34, ♦9000. /4 100. -HIT (I) 

102 CONTINUE 

DO 36 J=1,NST 
36 ALFA(J,I)=CS (JH-20, J) 

ALFA (JH-15,I)=-1.0 
GO TO 29 


K 


I 



35 IF (JH.EQ.23) GO TO 29 
BETA (I) =0.0 
DO 37 J=1,NST 
ALFA (J,I)=0.0 
DO 37 L=1,NST 


29 CONTINDE 


RETURN 


END 




I 


I 


I 


I 


I! 


I 
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SUBFOOTINE THAJ 

IMPLICIT REAl*8 (A - H, 0 - Z) 

BEAL*8 ID 

COMHON/CA/C(8,3,5) ,D(8,4,5) ,E(8) ,XI(3,5) , 01 (4 , 5) ,T I (8, 5) ,I0 (7) , 
1XD(3) 

COHMON/CB/A(7,7,20) ,BX(7,20) ,T (20) , VAL (80) , JM (80) , HIT (20) ,JHIT(20) 
COnnON/CO/ALPA(7,20) ,BETA(20) ,PAHnS(3,4) ,TS (20) ,XSTO(7,20) , 

1XLSTO (7,20) ,XLO (7) , IT (3) ,JM (5) ,NPHASE 
COMHON/CX/P(3,3.5) ,G{3,4,5) ,00(4,4,20) ,10(7,7) 

COMHON/NS/NSTATE 

COHMOM/TTPANS/PAR(20,20) , PA BIN (3,3) ,PAfiS (3,3) ,X(7) ,KPLAG,KHAX, 
1KOPT(20) ,KOUNT,NOPT (20) , NOPTS 
DIMENSION Y (7) ,H (7) 

KFLAG=0 

1 DO 2 I=1,NSTATE 

2 X(I)=X0(I) 

KOUNT=KOUNT+1 

KT=0 

DO 4 M=1,NPHASE 

84 CALL UPDATE(X,T (M) ,A(1,1,H) ,BX (1,M) ,NSTATE, 1) 

DO 104 I=1,NSTATE 

104 XSTO(I,M)=X (I) 

88 IP(JHIT (M) ,EQ. Q> GO TO 80 
S=BETA (M) 

DO 81 I=1,NSTATS 

81 S=S-ALFA (I,H) ♦X(I) 

IF(DABS (S/BETA (M)) .LE.. 0001) GO TO 83 
CALL VMULT(A(1,1,M) ,X,Y,NSTATB) 

SD=0.0 

DO 82 1=1,NSTATE 

82 SD=SD-ALFA (I,M) * (Y (I) ♦BX (I, M) ) 

T(M)^T (H)“S/SD 

KT=KT+1 

IP(KT.LE.10)GO TO 90 

KPLAG=1 

BET URN 

90 DTM=-S/SD 

CALL UPDATE(X,DTM,A (1,1, M) , BX (1,M) , NSTATE, 1) 

DO 105 I=1,NSTATE 

105 XSTO(I,M}-X(I) 

GO TO 88 

83 IF(H,GT.1)TS(M)=T(H)+TS(M-1) 

IF(M. EQ- 1) TS(1)=T(1) 

T(M+1)=TS(W + 1)-TS(M) 

80 CONTINUE 

IF (IT (1) .EQ.O) GO TO 4 
IF(K.EQ.0)GO TO 50 
DO 51 J=1,K 

CALL VMOLT(PHI,PAR (1,J) ,«,NSTATE) 

DO 51 1=1,3 
51 PAR(I, J)=W(I) 

50 DO 100 1=1,3 
100 PAR(I,M)=Y(I)+BX(I,M) 

DO 102 K=l,4 


IN=JM (K) 

IF (H. NE.IN) GO TO 102 
DO 85 1=1,3 
L=IT{I) 

IF (I.GT.IN) GO TO 85 

SDS=0.0 

KX=K-1 

DO 86 J=1,3 

SDS=SDS-ALFA (J , K) ♦ ( PA R( J ,L) -PAR{J,L + 1)) 

IF(KX.EQ.0)GO TO 86 
DO 103 JJ=1,KX 

IF(JK.LT.L.OR. JK.GE.10) GO TO 103 

SDS=SDS-ALFA (J,K)*PARMS (I, J J) ♦ (PAR ( J , JK) -P AR (.7, JK+ 1 )) 
103 CONTINUE ' 

86 CONTINUE 

PARMS (I,K)=-SDS/SD 
85 CONTINUE 
102 CONTINUE 
4 CONTINUE 

IF(IT (1) .EQ.O) GO TO 106 

87 DO 61 J=1,3 
K=IT(J) 

DO 61 1=1,3 

PARS (I,J)=PAR(I,K) -PAR(I,K+1) 

DO 61 M=1,4 
in=Jf1 (M) 

IF (IK. LT. K. OR.IH.GE.10) GO TO 61 

PARS(1, J) =PARS (I,J) +PARMS(J,M) ♦ (PAR (I, IK) -PAR (I,IH + 1) 1 
61 CONTINUE ^ f n 

106 CONTINUE 
RETURN 
END 
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SUBEOUTINE aPDATE( X,T, A , BX, ESTATE, L) 

IMPLICIT REAL*8 (A - H, 0 - Z) 

DIHENSION X(NSTATE) ,BX (ESTATE) , A (ESTATE, ESTATE) ,BS (10) ,BT(10) , 
1XS(10) ,XT(10) 

E=12 

IP(L.EQ.2) GO TO 7 
DO 6 1=1 , ESTATE 
BS(I)=BX (I) 

6 XS(I)=X(I) 

DO 2 K=1,E 

DO 3 1=1, ESTATE 

XT(I)=0.0 

BT(I)=0.0 

DO 3 J=1, ESTATE 

BT (I) =BT (I) *A (I, J) ♦T/DFLOAT (E-K + 2) *BS (J) 

3 XT(I)=XT (I) *A (I,J)*T/DPLOAT (E-K*1)*XS(J) 

DO 4 1=1, ESTATE 

BS (I) =BT (I) +BX(I) 

4 XS(I)=XT(I) ♦X(I) 

2 COETIEUE 

DO 5 1=1, ESTATE 

5 X(I)=XS(I)*BS(I)*T 
RETURE 

7 DO 1 1=1, ESTATE 
1 XS(I)=X(I) 

DO 8 K=1,E 
DO 9 1=1, ESTATE 
XT (I) =0.0 
DO 9 J=1, ESTATE 

9 XT(I)=XT(I)-A(J,I) ♦T/DFLOAT (N-K + 1) ♦XS (J) 

DO 10 1=1, ESTATE 

10 XS(I)=XT(I)+X(I) 

8 COETIEDE 

DO 11 1=1, ESTATE 

11 X(I)=XS(I) 

RETURE 

EED 
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